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1.0  INTRODUCTION 


The  Airframe  Propulsion  and  Weapons  Integration  Section  of  Wright  Laboratory  is  in  the 
process  of  expanding  its  computational  fluid  dynamics  (CFD)  capabilities.  The  charter  of  the 
Section  with  regard  to  CFD  and  propulsion  is  the  application  of  existing  CFD  tools  on  advanced 
inlet  and  nozzle  components  for  purposes  of  analysis,  evaluation,  and  design.  Recent  progress 
in  applications  on  various  nozzle  configurations  has  been  very  promising,  and  a  certain  level  of 
expertise  has  been  achieved.  However,  application  to  high  speed  inlet  components  has  not 
progressed  to  the  same  level.  Lack  of  experience  with  such  problems  has  led  to  limited  success 
and  is  the  primary  reason  for  slow  progress. 

An  extensive  literature  search  was  conducted  for  CFD  applications  to  high  speed  inlets 
and  inlet  related  flowfields.  The  range  of  materials  collected  include  two-dimensional  inviscid 
calculations  of  "unit"  problems  (compression  ramps,  shock  wave-boundary  layer  interactions,  etc.) 
to  three-dimensional  fully  turbulent  Navier-Stokes  analyses  of  real  inlet  test  hardware.  This 
information  dates  back  to  the  mid  1970's,  when  the  first  numerical  solutions  for  supersonic 
turbulent  flows  with  shock-boundary  layer  interactions  were  being  obtained.^'^  The  range  of 
numerical  applications  which  have  been  reviewed  is  shown  in  Fig  1-1.  Applications  to  purely 
axisymmetric  inlet  configurations  were  excluded  from  the  literature  search.  The  majority  of  the 
literature  covers  applications  in  the  Mach  range  2-7,  with  a  few  cases  for  Mach  number  greater 
than  12.  The  altitude  range  where  these  inlets  are  operating  as  part  of  an  airbreathing  propulsion 
system  is  generally  below  150,000  feet.  The  expected  flow  condition  for  this  regime  is  a  real 
gas  continuum.  Elevated  temperatures  in  stagnation  regions  and  boundary  layers  support  the 
possibility  of  vibrational  and  chemical  equilibrium  thermodynamics.  However,  finite  rate 
chemical  reactions  are  usually  not  considered  for  inlet  flows.  High  flight  Reynolds  number 
combined  with  shock  wave  induced  disturbances  guarantee  that  turbulent  boundary  layers 
predominate  the  viscous  interactions.  In  some  wind  tunnel  environments  where  Reynolds  number 
is  relatively  low  with  an  isolated  inlet  installation,  laminar  boundary  layers  may  extend  a 
significant  distance  aft  of  leading  edges.  Boundary  layer  transition  is  observed  to  occur  in  such 
cases. 

The  purpose  of  this  report  is  to  review  important  aspects  of  CFD  application  to  high  speed 
inlets  and  inlet  related  flowfields  based  on  the  data  in  the  literature  search.  Topics  which  are 
covered  include  governing  equations,  numerical  integration  algorithms,  grid  requirements, 
numerical  boundary  conditions,  and  turbulence  models.  When  possible,  the  results  of  numerical 
analyses  compared  to  experimental  data  are  shown  to  allow  an  assessment  of  the  CFD 
application. 
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Figure  1-1  Survey  Applications 
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2.0  GOVERNING  EQUATIONS 


The  application  of  CFD  techniques  to  any  aerodynamic  problem  starts  with  the 
identification  of  the  governing  equations  of  the  fluid  motion.  As  mentioned  in  the  introduction, 
high  speed  inlets  operate  in  the  compressible  flow  regime  usually  at  flight  Reynolds  numbers 
which  result  in  turbulent  flow.  The  most  general  and  applicable  equations  which  govern  this  flow 
are  the  time  dependent  Reynolds  averaged  compressible  Navier-Stokes  equations.  For  Cartesian 
velocity  components  (u,v,w),  the  Navier-Stokes  equations  are  written  in  strong  conservation  law 
form  as 


dU.  dE.  dF. 

dy  dz 


2.1 


The  conservation  variable  vector  is  U  =  {p ,  pu,  pv,  pw,  Ef.)  ^  ,  and  the  flux  vectors  are 


E  = 


Q  = 


pu 

pu2+p-T^ 

pv 

pu^-V 

pV^+p-T^ 

puv-V 

[E^.+P)  V-UXy^-VXyy-VrZy^+q^^ 

pw 

puw-x^ 
pvw-x^ 
pw^+p-x^^ 

(E^+p)  w-ux^-vx^-wz^^+gj 


2.2 


E^  =  p  e+ 


u^+v^+w^ 
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The  viscous  stress  terms  are 


-r  =  2  du_  0V_^\ 

x^=  I  (2|r-|H-|^, 


and  the  heat  flux  terms  are 


g,  =  -  (*+*,)  . 

=-(*+*,)  H, 

g,  =  . 


In  deriving  the  viscous  stresses,  it  is  assumed  that  the  turbulent  stresses  are  proportional  to  the 
mean  flow  strain  rates  in  the  same  manner  as  molecular  viscous  stresses.  This  is  referred  to  as 
the  Boussinesq  assumption,  and  is  discussed  in  more  detail  in  Section  2.2.  This  set  of  differential 
equations  are  of  course  based  on  the  conservation  principles  of  mass,  momentum,  and  energy. 
Their  integration  in  time  and/or  space  via  numerical  methods  is  the  basic  objective  of  CFD. 

In  order  to  make  the  Navier-Stokes  equations  tractable  for  numerical  integration,  they  are 
usually  non-dimensionalized  and  transformed  for  a  uniformly  spaced  computational  coordinate 
system.  The  transformation  is  pictured  graphically  in  Figure  2-1.  See  Reference  3  for  a 
complete  development  of  this  process.  The  transformation  simplifies  the  discretization  process 
and  allows  the  digital  computer  program  (or  "code")  to  systematically  store  and  retrieve  the 
necessary  information  quickly.  The  disadvantage  is  the  requirement  to  produce  a  highly 
structured  field  grid  which  can  be  a  time  consuming  interactive  endeavor. 

2.1  Gas  Properties 

In  addition  to  Equations  2.1  through  2.4,  the  gas  state  relationships  between  the 
thermodynamic  variables  (p,p,T,e,h)  are  needed  as  well  as  the  variation  of  the  molecular  transport 
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properties  (w,jt)  with  temperature.  The  complexity  of  these  relationships  depends  upon  the  gas 
temperature  range  of  the  flows  being  considered. 

At  the  pressures  and  temperatures  expected  in  inlet  flows,  air  behaves  according  to  the 

perfect  gas  law: 


p  =  pRT 


2.5 


where  R  is  the  gas  constant,  which  depends  upon  the  gas  composition. 

For  low  temperatures  (3K-800K),  air  is  also  considered  calorically  perfect,  where  the 
specific  heat  capacities  at  constant  volume  and  pressure  (c^,  and  c^)  are  constant.''  A  calorically 
perfect  gas  implies  all  of  the  following: 


e=c^T, 


h=CpT, 


R, 


Y-1' 


c=^ 
P  y-1 


2.6a 


The  values  of  c„,  c^,  R,  and  y  are  all  constant  under  the  calorically  perfect  gas  law. 

As  air  temperature  increases  above  about  800K,  diatomic  air  molecules  become 
vibrationally  excited.  In  this  case,  the  specific  heat  capacities  are  not  constant,  but  functions  of 
temperature: 

2.6b 

Cp  =  Cp{T) 


Hence,  e  and  h  are  related  to  temperature  through: 

de  =  CydT 

dh  =  CpdT  2.6c 

e  =  e(r) 
h  =  h{r) 

A  gas  is  considered  thermally  perfect  if  it  is  not  chemically  reacting  and  the  internal  energy  and 
enthalpy  are  functions  of  temperature  only.  The  gas  constant  R  remains  unchanged. 

The  molecular  viscosity  is  a  function  of  temperature  only  for  low  density  fluids  such 
as  air,  and  is  normally  determined  from  kinetic  molecular  theory  and  empiricism.^  One  relation 
is  the  well  known  Sutherland's  law, 

ji  /  2  7a 

jXq  "  [t,)  ~Ws 

where  the  subscripted  quantities  are  reference  values  and  S  is  a  constant.  Sutherland's  law  is 
accurate  for  air  to  within  2%  in  the  temperature  range  170K-1900K.  The  power  law. 


5 


n«0 .7 


2.7b 


JL  = 

V-o 


«)■■ 


is  slightly  less  accurate  (4%)  over  the  same  temperature  range/  Of  the  applications  investigated 
here  which  disclosed  the  particular  viscosity  relation  used,  all  but  one  used  Sutherland's  law. 
Agarwal,  Deese,  and  Peters*  is  the  one  exception  which  used  the  simpler  power  law.  No 
particular  reason  is  given  for  their  selection.  It  is  doubtful  that  one  could  Identify  any  difference 
in  inlet  CFD  solutions  obtained  with  these  different  molecular  viscosity  relations,  especially  since 
turbulence  dominates  the  viscous  phenomena  anyway.  The  fact  that  many  of  the  application 
references  do  not  even  disclose  which  relation  is  used  is  further  evidence  that  the  selection  is  not 
critical.  The  one  advantage  the  power  law  has  is  that  it  can  be  completely  non-dimensionalized; 
Sutherland's  law  requires  a  characteristic  temperature  value  for  non-dimensionalization. 

The  thermal  conductivity  k  is  also  a  function  of  temperature,  but  can  be  related  to 
molecular  viscosity  via  the  Prandtl  number: 

Pr  =  2.8 

k 


where  the  Prandtl  number  is  assumed  to  be  constant.  The  Prandtl  number  was  0.72  (air  at 
standard  conditions)  for  most  of  the  applications  reviewed  for  this  report.  A  value  of  0.73  was 
used  in  a  few  applications.  None  of  the  references  discuss  the  impact  of  choosing  one  value  over 
the  other;  the  result  would  be  on  the  order  of  one  percent  difference  in  thermal  conductivity. 
Again,  it  is  unlikely  that  any  difference  could  be  discerned  in  inlet  solutions. 

At  higher  Mach  numbers,  temperatures  at  stagnation  points  and  in  boundary  layers  may 
be  high  enough  to  cause  molecular  dissociation  and  other  chemical  reactions  (at  atmospheric 
pressure,  O2  dissociates  at  2000K  and  Nj  dissociates  at  4000K).'‘  These  are  energy  absorbing 
processes.  Under  these  conditions,  the  specific  heat  capacities  are  functions  of  temperature  and 
pressure  and  the  gas  chemical  composition  can  change.  A  high  temperature  gas  model  is  used 
for  the  gas  state  relations  and  the  transport  coefficients.  These  models  are  called  equilibrium  air 
models.  They  return  the  experimentally  observed  equilibrium  state  of  the  gas  as  a  function  of 
temperature  and  pressure.  The  molecular  transport  coefficients  are  also  determined  as  a  function 
of  temperature.  Equilibrium  air  models  are  generally  empirical  curve  fits  or  table  look-up 
routines,  and  as  such  require  more  computing  resources  than  the  simple  algebraic  expressions  of 
the  perfect  gas  law  and  calorically  perfect  gas  relationships.  None  of  the  applications  reviewed 
here  considered  non-equilibrium  chemistry.  It  is  usually  assumed  that  the  chemical  reactions  are 
fast  compared  to  the  inlet  flow  residence  time.  Finite  rate  chemistry  gas  models  do  exist, 
however,  and  are  often  used  for  hypersonic  nozzle  flows. 

References  7,  8,  and  9  have  studied  the  effects  of  different  values  of  y  in  perfect  gas 
calculations  as  well  as  full  equilibrium  air  models  on  the  NASA  P8  inlet  model. This  model, 
shown  in  Figure  2-2,  was  tested  at  Mach  7.4  and  continues  to  be  a  favorite  subject  of  CFD 
application.  The  freestream  conditions  were  air  at  Mach  7.4,  Reynolds  number  of  8.86  million 
per  meter,  and  total  temperature  of  81  IK.  The  model  is  a  mixed  compression  fixed  geometry 
inlet  for  scramjet  propulsion.  The  design  incorporates  a  single  external  ramp  shock,  and  a  cowl 
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physical  space 


computational  space 


Figure  2-1  Coordinate  Transformation  Schematic 


Figure  2-2  NASA  P8  Inlet  (Ref.  10) 


shock  which  reflects  weakly  off  the  ramp  shoulder.  No  boundary  layer  separation  was  observed 
at  this  interaction.  Instrumentation  was  primarily  to  map  to  the  flowfield,  rather  than  measure 
performance.  Wall  static  pressure  and  temperature  distributions  were  taken  on  the  ramp  and 
cowl.  Also  pitot  pressure  and  total  temperature  vertical  surveys  from  ramp  to  cowl  were  taken 
at  many  streamwise  stations  throughout  the  model.  Successfulness  of  CFD  applications  to  the 
P8  has  generally  been  based  upon  comparison  to  the  pitot  pressure  profiles. 

For  CFD  analysis,  Mach  7.4  is  high  enough  to  warrant  the  use  of  high  temperature  gas 
properties.  However,  the  freestream  total  temperature  of  the  wind  tunnel  was  well  below  flight 
conditions,  which  lead  the  authors  of  these  references  to  compare  numerical  results  using  perfect 
gas  versus  high  temperature  gas  models. 

Kunik,  Benson,  Ng,  and  Taylor^  compared  values  on  the  P8  inlet  for  y  of  1.4  and  1.38 
in  perfect  gas  calculations.  The  technique  of  using  y  <  1-4  with  a  perfect  gas  model  is  a 
convenient  way  to  investigate  first  order  effects  of  high  temperature  without  the  complexity  of 
the  full  equilibrium  air  model.  They  concluded  that  there  were  significant  differences  in  the  pitot 
pressures  obtained  and  in  the  location  of  the  shock  waves,  as  would  be  expected.  Their  results 
show  that  the  pitot  pressures  for  y  =  1.38  are  much  closer  to  the  experimental  data  as  compared 
to  the  y=1.4  resutls.  (The  results  are  not  shown  here  due  to  poor  quality  of  the  reference  copy.) 

Rhie  and  Stowers®  compared  results  on  the  P8  inlet  for  of  a  perfect  gas  model  with  y  = 
1.38  to  an  equilibrium  air  model.  They  concluded  that  the  equilibrium  air  model  did  not  show 
any  significant  difference  as  compared  to  the  perfect  gas  with  y  =  1.38.  They  did  not  perform 
calculations  with  y  =  1.4. 

Brenneis  and  Wanie’  performed  calculations  on  the  P8  inlet  with  perfect  gas  using  y  = 
1.4  and  1.38  and  also  an  equilibrium  air  model.  The  y  =  1.38  perfect  gas  calculation  matched 
the  experimental  pitot  pressure  profiles  best  overall  as  shown  in  Figure  2-3.  Note  that  all  three 
calculations  are  significantly  different  from  each  other,  indicating  that  the  pitot  pressure  is  very 
sensitive  to  changes  in  y.  In  contrast,  the  authors  mention  that  the  wall  static  pressure 
distributions  were  not  affected  by  the  gas  models.  The  authors  conclude  that  high  temperature 
gas  effects  are  present  in  their  solutions.  They  suggest  that  to  fully  compare  the  CFD  results  to 
test  data,  the  CFD  solution  would  have  to  include  the  flow  about  each  pitot  probe.  (In  dently, 
the  equilibrium  air  calculation  took  20%  more  CPU  time  than  the  perfect  gas  calculations.) 

Note  that  Brenneis  and  Wanie  observed  a  difference  in  their  results  between  y=1.38  and 
the  equilibrium  gas  model,  while  Rhie  and  Stowers  found  no  significant  difference.  The 
equilibrium  air  models  used  in  these  separate  references  are  from  different  sources,  which  might 
explain  the  conflicting  conclusions. 

It  is  obvious  from  the  above  discussion  that  high  temperature  gas  effects  are  first  order 
for  the  P8  inlet  test,  even  at  the  low  freestream  stagnation  temperature.  As  a  general  rule  of 
thumb,  high  temperature  gas  models  probably  should  be  considered  for  inlet  calculations  above 
a  Mach  number  of  about  6. 

2.2  Time  Averaged  Turbulent  Stress  and  Heat  Transfer 

The  turbulent  viscosity  (or  eddy  viscosity)  coefficient;/,  in  Equation  2.3  is  a  result  of  two 
important  developments  of  the  unsteady  Navier-Stokes  equations.  First,  it  is  assumed  that  flow 
variables  in  a  turbulent  flow  can  be  decomposed  into  a  time  averaged  mean  plus  a  fluctuating 


8 


component: 


<J)  =  4>+4>^  <1>  =  2.9 

At  tg 

where  (J)  could  be  any  of  the  flow  variables.  The  time  increment  At  is  large  compared  to  the 
fluctuations,  but  relatively  small  compared  to  the  time  constant  of  any  unsteadiness  in  the  mean 
flow.  This  is  depicted  in  Figure  2-4.  The  decomposed  flow  variables  are  substituted  into  the 
governing  equations  which  are  then  time  averaged.  In  this  process,  most  of  the  fluctuating  terms 
drop  out  by  definition  of  the  time  average  or  are  neglected.  The  extra  terms  that  remain  are  due 
to  correlations  of  the  fluctuating  velocity  components.  See  Reference  3  for  the  details  of  this 
derivation.  These  terms  are  referred  to  as  apparent  turbulent  stresses.  Second,  to  cast  the 
equations  in  the  form  given  in  Equations  2.1-2.4,  it  is  further  assumed  that  the  apparent  turbulent 
stresses  are  proportional  to  the  mean  flow  strain  rates  through  the  turbulent  eddy  viscosity 
coefficient.  This  is  called  the  Boussinesq  assumption  and  is  made  by  every  reference  dealing 
with  turbulent  flow  reviewed  herein. 

The  turbulent  coefficient  of  thermal  conductivity,  k,,  in  the  energy  equation  is  also  a  result 
of  the  same  developments  discussed  above.  The  apparent  turbulent  heat  flux  is  a  result  of 
correlations  of  the  fluctuating  velocity  components  and  fluctuating  temperature.  The  Boussinesq 
assumption  relates  the  apparent  turbulent  heat  flux  to  the  gradients  of  the  mean  temperature 
through  the  turbulent  thermal  conductivity  coefficient.  In  practice,  this  coefficient  is  related  to 
the  turbulent  viscosity  coefficient  by  a  turbulent  Prandtl  number: 

Prt  =  V^t^p/kt  2.10 

The  turbulent  Prandtl  number  Pq  is  taken  to  be  0.9  for  all  the  turbulent  applications  reviewed 
herein. 

There  are  two  different  commonly  used  time  averaging  procedures  for  the  governing 
equations.  The  first  is  a  strict  time  average  as  shown  in  Equation  2.9,  and  is  referred  to  as 
Reynolds  averaging.  The  other  is  a  mass-weighted  time  average,  or  Favre  average,  where  the 
time  mean  flow  variables  are  defined  as 


$  =  S  .  2.11 

P 

Some  authors  claim  advantages  to  mass-weighted  averaging.^’”  The  compressible  turbulent 
Navier-Stokes  equations  with  mass-weighted  variables  is  nearly  identical  to  the  incompressible 
turbulent  Navier-Stokes  equations.  According  to  Wilcox  and  Traci“,  this  is  consistent  with  a 
hypothesis  due  to  Morkovin  that  turbulence  in  equilibrium  behaves  in  a  locally  incompressible 
manner.  However,  flows  undergoing  rapid  expansions  or  compressions  are  not  in  equilibrium 
with  respect  to  turbulence.  Thus  compressibility  corrections  are  needed  in  the  turbulence  closure 
models.  Of  the  application  references  reviewed  here  dealing  with  turbulent  flows,  about  half  use 
mass-weighted  averaging  and  half  use  Reynolds  averaging.  No  direct  comparisons  were  made 
between  the  two  in  any  of  the  references. 


10 


t 


Figure  2-4  Unsteady  Mean  Flow  with  Fluctuating  Components 


•  Aspect  ratio  at  throat:  7.56 
at  exit:  2.80 

Figure  2-5a  Mixed  Compression  Inlet  Application  (Ref.  12) 
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2.3  Subsets  of  the  Navier-Stokes  Equations 

The  time  dependent  Navier-Stokes  equations  are  a  mixed  set  of  hyperbolic-parabolic 
partial  differential  equations^  which  govern  flow  behavior  for  initial-boundary  value  problems. 
They  are  solved  by  integrating  in  time  from  some  initial  condition  to  a  steady  state  (if  one  exists) 
with  influence  from  the  boundaries  of  the  domain  of  interest.  It  is  beneficial  to  determine  if 
appropriate  approximations  can  be  made  to  simplify  the  governing  equations  and  therefore  reduce 
computer  requirements  in  terms  of  memory  and  processing  time.  A  discussion  of  the  simplified 
equation  subsets  which  have  been  applied  to  high  speed  inlets  is  given  below,  including  thin  layer 
Navier-Stokes,  parabolized  Navier-Stokes,  and  Euler.  In  general,  all  of  the  problems  shown  in 
Figure  1-1  can  be  solved  using  state-of-the-art  Navier-Stokes  codes.  This  section  will  investigate 
which  flows  have  also  be  successfully  solved  with  the  simplifying  approximations. 

2.3.1  Thin  Layer  Navier-Stokes 

The  thin  layer  Navier-Stokes  equations  are  derived  by  neglecting  viscous  derivatives 
parallel  to  solid  walls.  The  usual  boundary  layer  order  of  magnitude  analysis  is  used  which 
retains  only  viscous  derivatives  in  a  direction  normal  to  the  wall.  This  can  be  accomplished 
easily  in  computational  coordinates  which  are  generally  oriented  with  the  body  surfaces.  The 
assumption  is  best  for  high  Reynolds  number  flows  where  the  boundary  layers  are  expected  to 
be  thin  relative  to  the  inviscid  core  flow  dimensions.  Unlike  the  boundary  layer  equations,  the 
thin  layer  equations  retain  all  other  terms  in  the  momentum  equation,  which  makes  them 
applicable  to  separated  and  reverse  flows  as  well.  Note  also  that  the  thin  layer  approximation 
neglects  both  molecular  viscous  stresses  and  turbulent  Reynolds  stresses  in  the  streamwise 
direction.®  The  advantage  over  the  full  Navier-Stokes  equations  is  a  reduction  in  the  number  of 
terms  needing  to  be  evaluated  by  the  computer  and  thus  reduced  CPU  time.  The  character  of  the 
equations  has  not  been  changed,  so  integration  in  time  is  still  required  for  a  solution.  The  thin 
layer  Navier-Stokes  equations  appear  to  be  attractive  for  high  speed  inlets  operating  in  a  high 
Reynolds  number  environment.  However,  one  might  question  the  ability  of  this  set  of  equations 
to  capture  the  complex  viscous  phenomena  of  shock  wave-boundary  layer  interactions,  shock 
induced  separations,  etc.,  which  are  typically  found  in  inlet  flows. 

Agarwal,  Deese,  and  Peters®  have  carefully  applied  a  thin  layer,  mass-averaged  code  to 
a  variety  of  shock  wave-boundary  layer  interactions,  including  2-D  and  3-D,  separated  and 
unseparated  flows.  The  primary  objective  of  their  work  was  to  study  the  role  of  turbulence 
closure  models  in  such  flows.  However,  they  show  through  order  of  magnitude  analysis  that  the 
thin  layer  approximation  is  appropriate  even  for  strong  interactions  (separated  flows)  when  the 
separation  bubble  is  "small."  In  general,  their  results  show  the  same  discrepancies  to 
experimental  data  as  full  Navier-Stokes  applications  of  the  same  flows.  Their  concluding  remarks 
indicate  that  the  thin  layer  approximation  is  not  affecting  the  results;  the  reason  for  discrepancies 
is  attributed  to  turbulence  models  (see  Section  6.0). 

Liou,  Hankey,  and  Mace^^  have  applied  a  thin  layer  Reynolds-averaged  code  to  a  2-D 
supersonic  mixed  compression  inlet  (Figure  2-5a).  The  computed  flow  exhibits  all  the  prominent 
features  observed  in  the  test,  including  shock  induced  separation,  and  the  natural  unsteadiness  of 
the  terminal  shock  in  the  diverging  channel.  Comparisons  to  experimental  data  show  the  same 
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reasonable  agreement  as  full  Navier-Stokes  calculations  of  other  mixed  compression  inlets  (Figure 
2-5b).  Again,  the  reason  for  discrepancies  is  not  attributed  to  the  thin  layer  approximation. 

Panaras  and  Stanewsky^^  have  applied  a  thin  layer  code  to  a  fin  generated  3-D  shock 
wave-boundary  layer  interaction.  No  discussion  is  given  for  the  suitability  of  the  thin  layer 
approximation  to  this  type  of  flow,  except  that  it  has  been  applied  to  other  supersonic  corner 
flows.  Their  results  capture  even  the  fine  details  of  the  flow  observed  in  the  corresponding 
experiment,  as  shown  in  Figures  2-6a  and  2-6b.  The  success  of  their  results  is  attributed  to  a 
turbulence  model  which  is  finely  tuned  to  this  type  of  flow. 

Unfortunately,  none  of  the  applications  show  a  direct  comparison  of  thin  layer 
approximation  to  full  Navier-Stokes  results.  Also,  the  authors  discussed  above  did  not  indicate 
how  much  computer  resources  were  saved  by  using  the  thin  layer  approximation.  It  appears  from 
the  above  discussion  that  the  thin  layer  approximation  can  be  applied  successfully  to  high  speed 
inlet  flows,  provided  the  flow  separations  are  not  too  severe. 

2.3.2  Parabolized  Navier-Stokes 

The  Parabolized  Navier-Stokes  (PNS)  equations  are  derived  by  neglecting  streamwise 
derivatives  in  the  viscous  terms  of  the  steady  Navier-Stokes  equations.  With  a  restriction  to 
supersonic  inviscid  core  flow,  the  governing  equations  reduce  to  a  set  which  is  parabolic  in  the 
streamwise  direction.  A  solution  to  the  PNS  equations  then  becomes  an  initial-boundary  value 
problem  which  is  integrated  by  marching  in  space  rather  than  time.  An  order  of  magnitude 
analysis^  shows  that  the  exclusion  of  the  streamwise  viscous  terms  only  becomes  appropriate  for 
hypersonic  flows  when  the  freestream  Mach  number  is  about  five  or  higher.  The  exclusion  of 
these  terms  eliminates  the  ability  to  compute  streamwise  separations,  although  crossflow 
separations  can  still  be  predicted  since  all  crossflow  terms  are  retained. 

Theoretically,  steady  hypersonic  flow  can  be  completely  solved  in  one  sweep  of  the 
streamwise  direction,  including  fully  coupled  inviscid  and  viscous  interactions.  Furthermore,  only 
the  solution  at  the  current  marching  step  is  stored  in  core  memory,  not  the  entire  flow  domain. 
This  has  obvious  major  savings  in  computer  processing  time.  For  these  reasons,  a  larger  number 
of  grid  points  can  be  used  to  resolve  the  boundary  layer  than  for  full  Navier-Stokes  calculations. 
To  start  a  PNS  solution,  the  flow  field  must  be  known  completely  at  some  initial  station  from 
which  the  streamwise  marching  begins.  This  can  be  the  freestream  for  configurations  with  sharp 
leading  edges,  or  a  Navier-Stokes  solution  for  a  blunt  leading  edge. 

Strictly  speaking,  the  streamwise  space  marching  integration  assumes  no  upstream 
influence  and  thus  requires  the  inviscid  flow  to  be  supersonic  with  no  negative  velocity 
components,  i.e.,  no  separation.  However,  any  streamwise  pressure  gradients  will  exhibit 
upstream  influence  through  the  subsonic  portion  of  the  boundary  layer.  This  causes  the  numerical 
problem  to  be  ill-posed  and  will  result  in  unstable  "departure  solutions"  unless  approximations 
are  made  to  the  pressure  gradient  term.  The  sublayer  approximation  assumes  that  the  streamwise 
pressure  gradient  is  nearly  constant  across  the  boundary  layer  in  high  Reynolds  number  flows, 
and  thus  can  be  evaluated  in  whole  or  in  part  in  the  supersonic  region  of  the  boundary  layer. 
Other  PNS  codes^'*  use  a  multi-pass  technique  which  uses  the  results  of  one  streamwise  sweep 
to  estimate  the  streamwise  pressure  gradient  for  the  subsequent  sweep.  The  sweeps  continue  until 
no  changes  are  seen.  This  process  is  sometimes  referred  to  as  global  pressure  relaxation.  Still, 


13 


Figure  2-5b  Thin-layer  N-S  Results  on  a  Mixed  Compression 
Inlet (Ref.  12) 
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Figure  2-6a  A  3-d  Shock-Boundary  Layer  Interaction  Application 
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Figure  2-6b  Results  of  Panaras  and  Stanewsky  (Ref.  13) 
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the  total  processing  time  for  a  PNS  solution  is  much  less  than  a  time  integrated  full  Navier- 
Stokes  solution. 

Applications  of  the  PNS  equations  to  high  speed  inlet  flows  are  numerous,  primarily 
because  of  the  comparative  ease  with  which  solutions  can  be  obtained.  Researchers  at  NASA 
Lewis  have  used  their  PNS  code  for  the  P8  inlet  in  two  and  three  dimensions,  3-D  glancing 
shock-boundary  layer  interactions,  and  other  mixed  compression  inlets  (limited  to  the  supersonic 
portion  only).  These  applications  are  discussed  below. 

Kunik,  Benson,  Ng,  and  Taylor’  comment  that  PNS  analysis  of  mixed  compression  inlets 
may  be  appropriate  because  viscous  phenomena  which  cannot  be  solved  by  PNS  (i.e.,  separation) 
would  lead  to  unstarts  and  poor  design  anyway.  They  conclude  that  it  is  necessary  from  a  design 
standpoint  to  predict  the  onset  of  separation,  but  not  the  details  of  separated  flow.  They  used  2-D 
PNS  analysis  to  study  effects  of  grid  spacing,  turbulence  models,  and  high  temperature  on  the 
P8  inlet.  Because  the  solutions  were  quickly  obtained,  an  optimization  of  these  parameters  could 
be  performed.  As  described  earlier,  they  were  able  to  obtain  good  agreement  with  the  pitot 
pressure  profiles.  A  3-D  PNS  calculation  was  also  performed,  but  a  good  comparison  to 
experiment  cannot  be  made  since  the  P8  test  model  was  not  sufficiently  instrumented  for  3-D 
features.  The  CFD  results  do  show  a  substantial  amount  of  three-dimensionality  created  by  the 
cowl  shock  interaction  with  the  sidewall  (not  shown  due  to  poor  quality  original). 

In  Reference  15,  an  AGARD  working  group  was  tasked  to  determine  the  current 
capabilities  of  CFD  applications  from  participating  US  and  European  agencies  on  several  inlet 
related  test  cases.  The  NASA  Lewis  researchers  applied  their  PNS  code  to  the  fin  generated 
glancing  shock-boundary  layer  interaction  test  case  and  the  P8  inlet  again.  The  results  were 
compared  with  test  data  and  a  number  of  full  Navier-Stokes  calculations  from  the  various 
agencies.  Very  little  details  of  the  applications  were  given;  the  focus  was  on  the  results  and 
comparisons.  The  PNS  code  used  in  Reference  15  was  the  same  as  the  one  used  by  Kunik  et 
al.’ 

For  the  3-D  glancing  shock-boundary  layer  interaction,  the  PNS  code  was  able  to  predict 
qualitatively  the  features  of  the  interaction.  Figure  2-7  shows  the  wall  pressure  comparisons  with 
experiment  and  several  other  Navier-Stokes  calculations.  The  PNS  code  exhibits  the  correct 
amount  of  upstream  influence,  as  shown  by  the  initial  pressure  rise  at  y=3.  However,  the 
pressure  in  the  plateau  region  is  over  predicted.  The  fact  that  the  PNS  code  can  predict  the 
upstream  influence  is  a  consequence  of  this  type  of  flow,  not  the  PNS  method.  The  ellipticity 
of  the  3-D  glancing  shock  wave-boundary  layer  interaction  acts  mostly  in  a  spanwise  direction, 
and  the  3-D  PNS  equations  retain  all  of  the  spanwise  terms  as  noted  previously.  Note  that  even 
the  Navier-Stokes  calculations  have  discrepancies  in  one  or  more  regions  of  the  wall  pressure 
plot.  The  pitot  pressure  profiles  and  a  yaw  angle  profiles  shown  in  Reference  15  illustrate  that 
the  PNS  code  does  as  well  as  the  full  Navier-Stokes  codes  for  these  parameters.  The  conclusion 
for  this  case  in  Reference  15  is  that  all  the  CFD  analyses  (including  NASA  Lewis'  PNS)  were 
applicable  to  this  type  of  flow,  but  that  user  involvement  in  modeling  decisions  (turbulence 
models,  boundary  conditions,  etc.)  will  have  a  strong  influence  on  the  results. 

For  the  P8  application  in  Reference  15,  ten  different  CFD  analyses  were  compared  against 
each  other  and  the  test  data.  The  plots  are  not  shown  here  because  the  number  of  curves  per  plot 
makes  them  difficult  to  read.  It  appears  that  the  PNS  results  in  this  reference  do  not  agree  with 
the  test  data  as  well  as  the  results  of  Kunik  et  al.’  Again,  no  details  of  the  PNS  solution  were 
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given  in  Reference  15.  Since  Kunik  et  af  conducted  a  careful  optimization  of  the  parameters 
affecting  the  solution,  it  is  concluded  that  their  results  represent  the  best  application  of  NASA 
Lewis'  PNS  code  to  the  P8  inlet. 

Kim,  Buggeln,  and  McDonald^*  have  applied  a  different  "parabolized"  Navier-Stokes 
method  to  a  2-D  shock  wave-boundary  layer  interaction  with  bleed  and  a  3-D  rectangular  inlet 
with  spillage.  The  method  involves  identification  of  a  primary  flow  direction,  aligning  a 
coordinate  with  that  direction,  and  neglecting  all  diffusion  along  that  coordinate.  The 
approximation  is  appropriate  for  high  Reynolds  number  flow.  The  subsonic  region  of  the 
boundary  layer  is  treated  with  additional  approximations  to  continuity  and  the  normal  momentum 
equation  based  on  the  reasoning  that  the  subsonic  region  will  be  thin.  This  has  resulted  in  a 
space  marching  method  which  does  not  exhibit  the  same  tendency  for  "departure  solutions"  as 
other  PNS  codes. 

Their  application  to  an  incident  shock  wave-boundary  layer  interaction  with  bleed  is 
shown  in  Figure  2-8.  For  this  case,  the  upstream  Mach  number  is  only  2.5,  which  would  usually 
be  considered  too  low  for  PNS.  Separation  is  avoided  in  this  interaction  by  a  boundary  layer 
bleed  region.  The  wall  pressure  comparisons  show  under  prediction  at  both  ends  of  the 
interaction,  which  the  authors  attribute  partially  to  lack  of  upstream  influence  in  the  space 
marching  procedure.  It  is  also  probable  that  the  bleed  boundary  conditions  have  some  effect,  per 
the  discussion  in  Section  4.2.  The  pitot  pressures  agree  qualitatively  with  the  experimental  data. 
The  authors  conclude  that  their  marching  procedure  is  applicable  to  this  type  of  flow  with  wall 
bleed.  However,  caution  is  advised  here  because  there  is  much  experimental  data  showing  that 
shock  induced  separation  and  its  control  by  bleed  is  sensitive  to  bleed  amount  and  location.  A 
situation  where  separation  is  not  controlled  (due  to  insufficient  bleed  rate  or  poor  location)  is  an 
eventuality  for  engineering  applications.  In  this  case  the  PNS  method  may  not  produce  a 
solution,  or  may  produce  an  erroneous  solution. 

Kim,  Buggeln,  and  McDonald^®  investigated  the  3-D  mixed  compression  inlet  to  determine 
PNS  applicability  to  inlets  with  significant  crossflow  caused  by  spillage  and  sidewall  interactions. 
No  test  data  was  available  for  this  case.  The  PNS  method  was  shown  to  be  capable  of  capturing 
the  complicated  3-D  spillage  flow  near  the  sideplate  lip  and  the  3-D  glancing  shock-boundary 
layer  interactions  along  the  sidewalls  caused  by  the  ramp  and  cowl  shocks. 

An  example  of  a  global  pressure  relaxation  technique  is  given  by  Yaghmaee  and 
Roberts.^"  Their  method  uses  forward  differencing  of  the  streamwise  pressure  gradient  to 
simulate  ellipticity  of  the  pressure  field  in  the  subsonic  boundary  layer  region,  with  multiple 
streamwise  sweeps  to  converge  on  a  solution.  This  method  is  also  sometimes  referred  to  as 
partially  parabolized  Navier-Stokes  (PPNS).  The  authors  have  applied  this  technique  to  subsonic 
flows,  terminal  normal  shocks,  and  shock  wave-boundary  layer  interactions  in  other  references. 
Note  that  their  implementation  allows  the  code  to  switch  from  single  pass  PNS  to  multi-pass 
PPNS  at  user  discretion.  This  is  convenient  since  one  can  march  along  with  PNS  until  a  strong 
elliptic  region  is  encountered  and  then  switch  to  PPNS.  In  this  reference,  they  have  applied  the 
technique  to  the  P8  inlet.  Some  of  their  results  are  shown  in  Figure  2-9.  They  have  about  the 
same  level  of  agreement  to  the  test  data  as  was  shown  in  Reference  15.  The  authors  concluded 
that  application  of  PPNS  in  the  shock-boundary  layer  interaction  region  did  not  make  a 
significant  difference  for  the  P8  because  the  interactions  were  relatively  weak. 
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Pr^ure  Distribution  on  the  Cowl 


Pressure  Distribution  on  the  Centerbody  Pressue  Profile  at  1.041  meters 


Figure  2-9 


P8  Application  Using  PPNS  (Ref.  14) 


2.3.3  Euler  Equations 

The  Euler  equations  are  the  inviscid  subset  of  the  unsteady  Navier-Stokes  equations.  They 
have  the  same  temporal  characteristics  as  the  Navier-Stokes  equations,  and  so  are  integrated  in 
time  to  yield  a  solution.  The  advantages  numerically  are  twofold.  First,  all  viscous  terms  are 
dropped  which  means  fewer  computer  operations  per  time  step.  Second,  resolution  of  thin 
boundary  layers  is  not  necessary.  Boundary  layer  resolution  involves  many  closely  spaced  grid 
points,  resulting  in  a  much  larger  number  of  grid  points  for  a  given  flow  domain  than  an  Euler 
solution.  The  disadvantage,  of  course,  is  that  important  viscous  phenomena  and  interactions  are 
not  modelled.  Shock  wave-boundary  layer  interactions  are  unavoidable  in  high  speed  inlets;  their 
control  is  a  crucial  design  issue  and  the  focus  of  continuing  research.  Thus  it  would  seem  that 
the  Euler  equations  have  limited  utility  in  high  speed  inlet  applications.  However,  the  state-of- 
the-art  design  procedure  for  inlets  is  to  apply  inviscid  analysis  to  a  conceptual  design  to 
determine  approximate  shock  wave  locations  and  strengths,  pressure  gradients,  etc.  Boundary 
layer  displacement  corrections  and  empirical  shock  wave-boundary  layer  interaction  techniques 
are  then  applied  to  modify  the  design  as  necessary.  The  Euler  equations  can  provide  a  quick, 
inexpensive  solution  for  design  purposes. 

There  are  two  examples  in  this  review  where  the  Euler  equations  have  been  used  for  inlet 
analysis/design.  Bissinger  and  Eberle^’  have  used  an  Euler  solver  to  study  the  flow  in  a  Mach 
5.65  mixed  compression  inlet,  as  shown  in  Figure  2-10.  From  a  design  standpoint,  they  have 
investigated  starting  enhancement  by  bleed,  Mach  number  perturbations,  and  3-D  shock  structure. 
The  Euler  analysis  was  capable  of  predicting  a  small  subsonic  pocket  near  the  ramp  shoulder 
which  would  no  doubt  cause  an  unstart  in  a  viscous  flow.  Also,  the  3-D  solution  shows  the 
effects  of  the  sidewalls  and  the  resulting  curvature  of  the  ramp  shocks.  This  type  of  information 
is  needed  to  modify  the  design  to  avoid  unstart  and  reduce  sidewall  spillage. 

Fujimoto,  Niwa,  and  Sawada^®  have  obtained  an  Euler  solution  on  a  Mach  2.5  mixed 
compression  inlet  to  verify  their  method-of-characteristics  design  process.  The  results  shown  in 
Figure  2-11  indicate  the  Euler  and  method-of-characteristics  are  closely  matched.  However,  a 
Navier-Stokes  solution  of  the  same  configuration  shows  that  shock  induced  boundary  layer 
separation  will  cause  unstart  unless  substantial  bleed  is  used.  This  result  illustrates  the  usual 
limitation  of  Euler  analysis  in  inlet  design,  where  viscous  effects  can  seriously  alter  the  expected 
flow  patterns. 

Another  area  where  the  Euler  equations  are  useful  is  in  the  analysis  of  terminal  shock 
stability  and  unsteady  operational  phenomena.  Many  references  (17,18,19,20)  attempting  to  solve 
mixed  compression  inlets  with  terminal  normal  shocks  indicate  that  shock  position  is  highly 
sensitive  to  subsonic  exit  boundary  pressure,  proximity  of  the  exit  boundary  to  the  shock,  and 
flow  solver  numerical  integration  parameters.  These  are  numerical  modeling  issues  not  directly 
related  to  viscous  phenenoma  and  can  be  efficiently  investigated  with  Euler  analysis.  Unsteady 
phenomena  being  investigated  are  inlet  unstart,  buzz,  and  hammershock.  These  events  can  be 
triggered  by  shock  induced  boundary  layer  separation,  freestream  disturbances,  or  engine  surge. 
Shock  induced  separation  is  of  course  a  viscous  phenonmenon,  but  freestream  disturbances  and 
engine  surge  can  occur  independently  of  the  inlet.  The  Euler  equations  are  applicable  because 
once  the  event  is  initiated,  the  inlet  flow  is  dominated  by  inviscid  transient  shock  wave  motion. 
Unstart,  buzz,  and  hammershock  events  have  been  known  to  cause  structural  damage,  and 
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estimates  of  potential  inlet  duct  over-pressures  are  needed  for  structural  design  purposes. 

Hindash,  Bush,  and  Cosner^’  applied  the  Euler  equations  to  the  classic  1-D  shock  tube 
problem  to  verify  their  accuracy  for  analysis  of  hammershock  events.  They  found  that  the  shock 
wave  speed  and  strength  (pressure  rise)  could  be  computed  within  5%  of  the  analytical  solution. 
Further  analysis  of  a  fighter  inlet  duct  hammershock  event  is  shown  in  Figure  2-12.  The  pressure 
rise  along  the  duct  is  compared  to  flight  test  data  where  the  maximum  pressure  at  several  duct 
locations  was  recorded  as  the  shock  passed  forward.  The  computational  results  were  found  to 
be  within  5%  of  the  test  data.  The  authors  conclude  that  this  is  sufficient  for  structural  design 
purposes. 
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Euler  Analysis  of  Inlet  Duct  Hammershock  Event 
(Ref.  19) 
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3.0  NUMERICAL  INTEGRATION  SCHEME 


The  focus  of  scheme  development  is  the  discretization  and  linearization  of  the  nonlinear 
governing  equations  over  small  increments  of  time  and  space,  and  the  numerical  procedures  used 
to  solve  the  resulting  algebraic  system.  The  objective  is  to  obtain  a  scheme  which  solves  the 
system  of  equations  in  an  efficient  and  stable  manner.  This  section  is  meant  to  give  a  brief 
general  assessment  of  the  state-of-the  art  from  an  application  viewpoint  and  summarize  the 
information  obtained  from  the  literature  survey.  A  complete  description  of  all  the  techniques 
being  used  is  beyond  the  scope  of  this  report. 

The  critical  issue  of  numerical  schemes  for  supersonic  flow  is  the  accuracy  with  which 
shock  waves  are  resolved.  Note  that  inlet  applications  are  not  unique  in  this  regard;  schemes  that 
work  for  general  supersonic  problems  will  work  for  inlet  applications.  The  classical  numerical 
methods  use  finite  difference  representations  of  the  governing  equations.  Because  the  difference 
equations  are  derived  from  Taylor  series,  their  solutions  must  be  continuous.  Thus  shock  waves 
are  represented  by  steep  gradients,  rather  than  true  discontinuities.  Typical  results  for  older  codes 
using  classical  first  and  second  order  differencing  will  have  dissipative  smearing  and  dispersive 
oscillations  near  steep  gradients  due  to  the  truncation  errors  of  the  Taylor  series  approximations. 
Artificial  dissipation  or  damping  is  usually  necessary  to  control  the  dispersive  oscillations  for 
stability.  Flux  splitting  techniques  provide  directional  differencing  that  represents  the  natural 
wave  propagation  directions  of  the  governing  equations.  This  helps  reduce  dissipative  errors  and 
steepen  the  shock  gradients.  Two  classical  methods  have  been  predominant  in  the  CFD  codes 
applied  to  inlets  according  to  the  literature  search.  MacCormack's  predictor-corrector  method  in 
various  explicit  and  implicit  formulations  has  been  used  for  unit  problems  and  complete  3-D  inlet 
calculations.  The  Beam- Warming  implicit  method  is  implemented  in  the  familiar  PARC  series 
of  codes,  and  most  major  airframe  contractors  have  proprietary  versions.  These  codes  are 
routinely  used  for  inlet  applications  of  all  types. 

The  most  recent  scheme  developments  are  able  to  blend  low  order  methods  near 
discontinuities  with  high  order  methods  in  smooth  regions  of  the  flow.^^  These  developments  are 
based  on  finite  volume  formulations  of  the  governing  equations.  All  of  the  following  methods 
were  encountered  in  the  most  recent  inlet  application  references:  various  upwind  schemes, 
Riemann  solvers,  total  variation  diminishing  (TVD)  schemes,  and  monotone  upstream  centered 
schemes  (MUSCL).  These  advanced  schemes  have  proven  to  be  fully  capable  of  accurately 
predicting  the  inviscid  shock  structures  of  any  practical  inlet  geometry.  Scheme  development  has 
not  been  a  restrictive  technology  issue  for  inlet  applications  for  several  years. 

Efficiency  improvements  will  always  be  in  demand  as  GFD  is  applied  more  for  design. 
Minimizing  the  number  of  iterations  (time  or  pseudo-time  integration  steps)  to  converge  to  a 
steady-state  condition  is  critical  in  determining  the  number  of  parameters  that  can  be  investigated 
in  a  design  application.  Again,  inlet  applications  are  not  unique  in  this  regard.  The  latest 
methods  which  have  been  applied  are  localized  time  stepping,  Newton  iteration,  Runge-Kutta 
iteration,  and  multi-grid  techniques. 
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4.0  BOUNDARY  CONDITIONS 


Since  the  Navier-Stokes  equations  and  the  subsets  are  initial-boundary  value  problems, 
the  boundary  conditions  define  their  unique  solutions.  Although  mathematical  analysis  indicates 
the  number  and  type  of  boundary  conditions  required  for  various  physical  boundaries,  there  is  no 
rigorous  theory  for  numerical  implementation  of  boundary  conditions.  Heuristic  analysis  is  often 
used  to  develop  and  test  numerical  techniques  which  exhibit  the  appropriate  physical  behavior. 
Equations  with  hyperbolic  or  elliptic  characteristics  are  particularly  sensitive  to  boundary 
conditions  because  disturbances  and  errors  are  propagated  to  other  regions  of  the  domain.^ 

Boundary  conditions  for  high  speed  inlet  flows  are  more  complicated  than  most  other 
applications.  Inlet  flows  have  multiple  supersonic  and  subsonic  regions  in  the  same  domain,  solid 
walls  with  various  thermal  constraints,  bleed  and/or  blowing  regions,  external  spillage,  etc. 
Treatment  of  boundary  conditions  for  inlet  applications  is  extremely  important  to  yield 
meaningful  results.  A  discussion  of  boundary  condition  types  encountered  in  this  survey  is  given 
below. 

4.1  Freestream/Inflow  Conditions 

This  is  the  easiest  condition  to  apply  for  high  speed  inlets  and  related  flows.  Since  the 
flow  is  supersonic,  characteristic  theory  dictates  that  the  freestream  boundary  should  not  be 
affected  by  other  points  in  the  field.  Thus  it  is  appropriate  to  fix  the  flow  variables  at  the 
freestream  boundary  in  time  and  space.  The  only  requirement  is  that  the  specified  conditions 
fully  define  the  conservation  variable  vector  U  at  every  point  on  the  boundary.  The  condition 
can  be  specified  as  total  conditions  or  static  conditions  or  combinations  thereof  which  can  be 
manipulated  with  the  usual  gas  dynamic  relations. 

Most  applications  to  actual  inlet  geometries  have  uniform  freestream  conditions. 
Specification  of  a  freestream  Mach  number,  total  pressure,  and  total  temperature  are  sufficient 
to  define  the  state  of  the  gas  on  the  entire  boundary.  Some  applications  specify  upstream 
boundary  flow  profiles  in  order  to  avoid  computing  the  flow  outside  the  region  of  interest.  Even 
the  subsonic  portion  of  the  boundary  layer  is  elliptic,  the  freestream  boundary  can  be  placed  far 
enough  upstream  that  the  pressure  gradient  will  not  be  influenced  by  events  downstream.  Then 
the  upstream  profile  can  be  held  fixed  without  violating  the  physics  of  the  problem.  For 
example,  shock  wave-boundary  layer  interactions  have  a  limited  region  of  upstream  influence. 
Because  the  flow  ahead  of  this  region  is  not  affected  by  the  interaction,  it  is  not  necessary  to 
solve  for  the  developing  boundary  layer  far  upstream  of  the  interaction.  Specification  of  the 
profile  has  been  accomplished  by  three  different  methods  in  this  survey:  1)  computing  a  flat 
plate  boundary  layer  prior  to  computing  the  solution  on  the  region  of  interest,  2)  assuming  a 
turbulent  profile  anchored  with  limited  test  information  (skin  friction,  displacement  thickness, 
etc.),  or  3)  matching  a  measured  pitot  pressure  profile.  References  2,  13,  16,  22,  23,  24,  25,  26, 
and  27  have  all  used  these  techniques  successfully  for  2-D  and  3-D  interaction  problems. 

For  turbulent  calculations  with  two-equation  turbulence  models,  the  turbulent  conservation 
variables  also  need  to  be  specified  at  the  inflow  boundary.  Description  of  the  models  and  details 
of  the  boundary  conditions  that  have  been  used  is  given  in  Section  6.2. 
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4.2  Wall  Boundaries 


Solid  walls  have  velocity  conditions  and  heat  transfer  conditions  in  viscous  flow.  The 
velocity  condition  is  the  no-slip  condition,  i.e.  the  fluid  velocity  at  the  wall  is  the  same  as  the 
wall  velocity.  None  of  the  applications  reviewed  here  had  moving  walls,  so  the  velocity  at  the 
wall  is  simply  set  to  zero.  The  heat  transfer  condition  can  be  adiabatic,  a  specified  heat  transfer 
rate,  or  a  specified  temperature.  If  the  heat  transfer  is  specified  or  adiabatic,  then  the  wall 
temperature  is  determined  as  part  of  the  solution,  and  vice  versa.  Numerically  there  is  no 
difficulty  with  either  condition,  and  many  codes  can  handle  both.  A  specified  wall  temperature 
usually  remains  fixed  in  time  but  can  vary  spatially.  A  specified  heat  transfer  rate  involves  a 
simple  difference  approximation  to  the  temperature  gradient  normal  to  the  wall. 

Although  the  heat  transfer  boundary  conditions  are  easy  to  implement  numerically,  it  is 
challenging  to  specify  conditions  that  realistically  model  a  particular  physical  situation.  A 
realistic  condition  depends  upon  the  heat  transfer  of  the  inlet  structure  which  can  be  non-uniform 
and  time  dependent.  For  test  models  in  high  Mach  blow-down  wind  tunnels,  a  constant  cold  wall 
temperature  is  a  good  approximation  because  the  model  may  not  heat  up  appreciably  over  the 
short  duration  of  the  run.  The  P8  inlet  is  a  good  example  here.  From  Reference  10,  the  typical 
run  time  is  only  1-4  minutes  and  the  model  leading  edges  were  water  cooled.  The  wall 
temperature  increased  from  the  ambient  condition  by  15%,  at  most,  over  a  20  second  time  span. 
Brenneis  and  Wanie’  performed  calculations  on  the  P8  with  constant  temperature  and  adiabatic 
wall  boundary  conditions.  The  results  for  the  adiabatic  case  were  not  in  good  agreement  with 
the  test  data  and  the  solution  predicted  a  small  separation  bubble  not  observed  in  the  experiment. 
References  28  and  29  also  used  constant  wall  temperature  conditions  for  the  P8  based  on  the 
reported  test  data.  For  continuous  flow  tunnels  or  flight  conditions,  the  steady-state  wall 
temperatures  and/or  heat  transfer  rates  are  rarely  known  in  sufficient  detail  for  a  full  specification 
of  all  wall  boundaries.  From  the  reference  survey,  the  rule-of-thumb  is  to  use  adiabatic 
conditions  for  test  models  in  continuous  flow  tunnels.  This  is  probably  a  reasonable  estimate  for 
uncooled  test  models  since  heat  conduction  away  from  the  model  is  limited.  Also,  for  geometries 
which  have  no  associated  wind  tunnel  test  or  heat  transfer  analysis,  the  adiabatic  condition  has 
been  used  most  frequently.  The  adiabatic  condition  at  least  yields  the  upper  limit  of  wall 
temperature,  which  might  indicate  if  cooling  is  needed  for  a  piece  of  hardware. 

Wall  pressure  determination  requires  consideration  of  the  wall  normal  momentum 
equation.  One  condition  is  to  set  this  gradient  to  zero. 


if  I 
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=  0 
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per  the  usual  thin  boundary  layer  assumption  for  high  Reynolds  number  flows.  This  is  the  most 
common  condition  from  the  survey  results  and  it  is  easy  to  implement:  the  pressure  at  a  wall 
point  is  set  equal  to  the  pressure  at  an  adjacent  point  normal  to  the  wall.  This  is  may  not  be  a 
good  approximation  for  hypersonic  flows  because  the  thin  boundary  layer  assumption  loses 
validity Large  magnitude  local  wall  curvature  also  tends  to  produce  a  non-zero  normal  pressure 
gradient. 

Alternatively,  the  steady  state  normal  momentum  equation  written  explicitly  for  the 
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pressure  gradient  is  used.  Since  the  wall  velocity  components  are  zero  (assuming  no  bleed  is 
present),  all  convective  terms  drop  out  and  the  normal  momentum  becomes: 


dP\ 

dr\ 


viscous  stress  derivatives 


4.2 


The  right-hand  side  in  complete  form  would  contain  all  the  viscous  stress  terms  in  Equation  2.3, 
transformed  into  a  wall  aligned  coordinate  system.  Some  of  the  reference  applications  neglect 
viscous  derivatives  in  the  tangential  directions  in  favor  of  the  normal  direction,  which  results  in 
a  reduced  normal  momentum  equation.  These  simplifications  are  tending  to  the  limiting  case  of 
Equation  4.1.  Theoretically,  a  reduced  normal  momentum  equation  should  be  more  accurate  than 
the  zero  normal  pressure  gradient  assumption  since  more  terms  are  retained.  Shang  et  al.^°’^^  have 
used  both  methods  for  shock  wave-boundary  layer  interactions  and  supersonic  compressing  comer 
flows  with  practically  identical  results.  These  results  were  for  solid  walls  with  no  bleed.  When 
wall  transpiration  is  involved,  the  form  of  the  wall  normal  momentum  equation  becomes  more 
important,  as  discussed  in  the  next  section. 

Other  methods  are  also  being  used  for  the  wall  pressure  gradient.  In  the  sublayer 
approximation  for  the  PNS  calculations  of  Kim  et  al.^®,  a  centrifugal  force  balance  equation  is 
used: 


dP\  ^  4.3 

di\  R 

where  V,  is  the  tangential  velocity  at  the  grid  point  adjacent  to  the  wall  and  R  is  the  local  radius 
of  curvature.  Liou  et  al.^^  use  an  unsteady  characteristics  method  to  determine  wall  pressure 
directly  in  the  inviscid  operator  of  their  numerical  integration  scheme. 

For  Euler  solutions,  wall  boundary  conditions  require  that  the  velocity  vector  be  tangent 
to  the  wall  since  mass  flux  can  not  pass  through  solid  walls.  This  is  more  tedious  to  implement 
numerically  than  the  zero  velocity  condition  for  viscous  flows.  None  of  the  Euler  application 
references  discuss  in  detail  how  wall  tangency  is  implemented.  Anderson,  Tannehill,  and 
Fletcher^  outline  two  methods  in  their  text  book:  a  reflection  procedure  which  involves  mirroring 
the  normal  velocity  component  and  thermodynamic  variables  about  the  wall  boundary  to  extra 
grid  points  beyond  the  boundary,  and  Abbett's  method  which  introduces  isentropic  compression 
or  expansion  waves  to  turn  a  misaligned  velocity  vector  parallel  to  the  wall. 

Accurate  specification  of  wall  normal  gradients  for  pressure  or  temperature  can  be  affected 
by  the  grid  topology.  The  normal  gradients  evaluated  as  one-sided  difference  approximations 
imply  that  grid  lines  in  the  r|-direction  are  normal  to  the  wall  everywhere,  a  subtle  detail  not 
often  discussed  in  application  references.  Grid  orthogonality  is  not  required  to  obtain  a  solution, 
but  too  much  deviation  from  orthogonality  will  produce  unacceptable  truncation  errors.^^^^ 
Thompson^^  indicates  that,  in  general,  near  orthogonality  at  the  boundaries  should  be  achieved. 
None  of  the  application  references  reviewed  here  investigated  the  effect  of  grid  orthogonality. 
Other  grid  considerations  are  discussed  in  Section  5.0. 
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4.3  Bleed  Boundary  Conditions 

Boundary  layer  bleed  subsystems  are  used  to  control  shock  boundary  layer  interactions 
in  high  speed  inlets.  Porous  wall  plates  and  slots  function  to  remove  all  or  part  of  the  low 
energy  air  near  the  wall  by  ducting  it  to  a  low  pressure  plenum  and  discharging  it  overboard. 
Bleed  location  and  amount  are  critical  design  parameters  for  control  of  shock  induced  separations 
and  overall  performance  of  the  inlet.  Bleed  boundary  conditions  are  therefore  an  important  aspect 
of  computational  analysis  of  realistic  inlet  systems. 

For  porous  wall  bleed,  it  has  been  impractical  to  model  the  individual  bleed  ports  in  a 
CFD  solution.  The  large  number  and  small  scale  of  the  bleed  ports  would  result  in  prohibitively 
large  grids  and  CPU  time.  Instead,  a  mass  flux  is  determined  which  is  allowed  to  transpire 
through  the  wall  at  the  discreet  grid  points  in  the  bleed  region.  The  most  common  approach  is 
to  specify  a  uniform  mass  flux  distribution  across  the  bleed  region  based  on  total  mass  flow  rate 
at  the  bleed  plenum  exhaust.^'^^’  This  is  straightforward  for  computing  solutions  on  experimental 
hardware  where  the  total  bleed  mass  flow  rate  is  measured.  For  design  purposes,  the  constant 
flux  specification  is  based  on  a  desired  amount  of  boundary  layer  removal.  The  amount  is 
specified  as  a  percentage  of  the  mass  flow  in  the  boundary  layer  or  a  percentage  of  the  total  inlet 
mass  flow  rate.  Trial  and  error  procedures  are  necessary  since  details  of  the  boundary  layer  and 
inlet  flow  are  not  known  a  priori.  Rizetta^^  experienced  difficulty  in  obtaining  solutions  when 
the  wall  velocity  changes  abruptly  from  zero  to  the  corresponding  bleed  normal  velocity  between 
adjacent  grid  points.  A  ramping  function  has  been  used  to  gradually  increase  the  mass  flux  from 
zero  to  the  uniform  level  over  a  few  grid  points.^®’^^ 

The  uniform  mass  flux  bleed  boundary  is  convenient  when  the  only  information  available 
is  the  total  bleed  mass  removed.  It  does  not  exhibit  the  true  nature  of  the  flow  in  the  bleed 
region  because  the  mass  flux  is  generally  not  uniform,  especially  when  shocks  fall  across  the 
region.  Experimental  results  indicate  that  the  bleed  mass  flux  increases  with  wall  pressure.^  The 
plenum  beneath  the  porous  plate  will  have  elliptic  upstream  influence,  and  strong  pressure 
gradients  across  the  plate  can  result  in  reduced  bleed  flux  or  even  recirculation  in  the  forward 
parts  of  the  bleed  region.  Some  applications  have  specified  a  non-uniform  bleed  flux  distribution 
based  on  more  detailed  experimental  information.^®  '*®"'*^  This  is  an  improvement  only  if  the  final 
wall  pressure  distribution  predicted  by  the  CFD  solution  matches  that  in  the  experiment.  If  the 
wall  pressure  is  different  than  the  experiment,  then  the  mass  flux  distribution  is  not  consistent 
with  the  experimental  conditions. 

There  are  current  ongoing  efforts  to  compute  the  small  scale  flowfields  in  isolated  bleed 
ports.  The  hope  is  that  by  understanding  the  detailed  aerodynamics  of  single  ports,  more  realistic 
bleed  flow  models  can  be  developed  to  improve  boundary  conditions  for  inlet  calculations  on  the 
macroscopic  level.  Hamed  et  al.'*^  '*^  have  computed  an  isolated  bleed  slot  on  a  flat  plate  with  and 
incident  shock  wave.  Their  results  (Figure  4-1)  show  that  the  flow  in  the  slot  acts  like  a 
convergent-divergent  passage,  which  becomes  choked  when  the  plenum  pressure  is  reduced  below 
a  certain  value.  The  velocity  distributions  show  that  the  normal  velocity  is  loaded  towards  the 
trailing  edge  of  the  slot  opening,  behind  the  incident  shock.  Also  the  tangential  velocity 
component  is  not  zero  (as  is  normally  assumed)  and  actually  reverses  to  an  upstream  direction 
near  the  trailing  edge. 

In  all  but  two  of  the  applications  reviewed  here,  the  wall  tangential  velocity  is  simply 
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specified  to  be  zero  for  lack  of  additional  information.  In  one  exception,  Hunter  et  al'*^  have  used 
additional  equations  assuming  constant  entropy  along  streamlines  in  the  bleed  region  which 
allows  the  wall  density  and  tangential  velocity  to  evolve  as  part  of  the  solution.  The  normal 
velocity  is  still  determined  from  a  specified  mass  flux.  Shigematsu  et  al'*^  have  used  this  same 
approach  for  a  Mach  3  mixed  compression  inlet.  None  of  the  applications  references  discuss  the 
impact  of  the  bleed  tangential  velocity  component  on  the  accuracy  of  results  obtained. 

Researchers  have  also  been  attempting  to  model  porous  wall  bleed  more  realistically  by 
coupling  the  bleed  mass  flux  with  the  static  pressure  in  a  fictitious  plenum,  without  actually 
computing  the  plenum  flow.  The  basic  method  is  to  use  the  pressure  difference  between  wall 
and  the  plenum  to  determine  the  mass  flux  through  the  bleed  openings.  Losses  are  incurred  via 
empirically  derived  flow  coefficients  (discharge,  velocity,  total  pressure  loss,  etc.)  for  the  bleed 
ports.  Abrahamson  and  Brower‘S  developed  an  initial  approach  by  treating  the  bleed  ports 
analytically  as  isentropic  convergent  passages  which  may  be  choked  or  unchoked.  The  wall 
tangential  velocity  is  assumed  to  be  zero  but  the  normal  velocity  is  determined  from  the  local 
wall  pressure,  a  specified  plenum  pressure,  and  an  empirical  velocity  coefficient.  The  velocity 
coefficient  at  each  grid  point  was  adjusted  to  match  the  experimentally  measured  mass  flux 
distribution  in  an  incident  shock-boundary  layer  interaction  test  case.  The  end  result  is  still  a 
specification  of  a  measured  mass  flux  distribution  by  an  indirect  method.  They  admit  that  this 
is  not  a  complete  verification  of  the  method  and  that  the  velocity  coefficient  will  be  a  function 
of  Mach  number,  Reynolds  number,  port  geometry,  etc.  This  method  does  have  good  potential 
assuming  empirical  relations  can  be  developed  for  the  bleed  port  flow  coefficients  across  a  wide 
range  of  conditions. 

Benhachmi,  Greber,  and  Hingst^  conducted  a  combined  experimental  and  numerical 
investigation  of  a  porous  flat  plate  with  an  incident  2-D  shock  wave.  The  porous  area  was  a 
woven  mesh  material  backed  by  a  honeycomb  as  opposed  to  a  solid  plate  with  holes.  Their 
experiments  included  measurements  of  surface  pressure  and  temperature  across  the  bleed  region 
as  well  as  mass  flux  distributions  entering  the  plenum.  The  CFD  analysis  included  an  empirical 
bleed  method  relating  the  mass  flux  distribution  to  the  wall  static  pressure  and  a  specified  plenum 
pressure.  The  air  resistance  to  the  porous  mesh  was  calibrated  prior  to  the  incident  shock  test 
conditions  to  yield  an  equation  of  the  form: 
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where  Re,  is  the  Reynolds  number  based  on  the  mesh  thickness,  and  A  and  B  are  empirical 
constants.  It  is  assumed  though  not  completely  clear  that  V  is  the  normal  wall  velocity 
component  resulting  from  the  pressure  difference  AP.  The  tangential  velocity  is  still  taken  to  be 
zero  in  the  CFD  analysis.  Equation  4.4  was  used  as  a  boundary  condition  for  cases  involving 
incident  shocks  created  by  wedge  angles  of  four,  six,  and  eight  degrees.  A  typical  result  of  the 
mass  flux  and  wall  pressure  is  shown  in  Figure  4-2  for  the  six  and  eight  degree  cases.  The 
perturbations  in  the  measured  mass  flux  in  the  bleed  region  are  due  to  roughness  of  the  mesh 
material,  small  steps  in  the  transition  from  non-porous  to  porous  walls,  and  possible  resonance 
in  the  plenum.  The  CFD  results  have  smooth  oscillations  probably  due  to  insufficient  streamwise 
grid  resolution  according  to  the  authors;  these  should  not  be  construed  as  the  same  perturbations 


31 


in  the  test  data  created  by  the  mesh  roughness.  The  authors  comment  that  bleed  region  roughness 
and  pressure  waves  emanating  from  the  plenum  are  impractical  to  model  in  a  discretized 
numerical  analysis.  Qualitatively,  the  calculated  bleed  flux  increases  with  the  wall  pressure, 
which  was  the  desired  result  of  this  bleed  boundary  condition.  The  predicted  bleed  mass  flows 
using  the  calibrated  loss  coefficients  are  within  seven  to  twelve  percent  of  the  measured  flow 
rates,  which  the  authors  suggest  is  good  for  engineering  applications  considering  the  simplicity 
of  the  model. 

The  bleed  boundary  conditions  are  also  dependent  on  the  wall  pressure  gradient 
assumptions  discussed  in  Section  4.2.  This  is  due  to  the  fact  that  the  streamlines  in  the  boundary 
layer  are  turned  toward  the  wall  when  suction  is  applied,  and  streamline  curvature  cannot  occur 
without  a  pressure  gradient.  Thus  assuming  zero  normal  pressure  gradient  is  not  physically 
correct  for  bleed  boundaries,  even  though  it  is  used  quite  often  for  simplicity.  Abrahamson  and 
Brower"**  applied  both  the  zero  normal  pressure  gradient  condition  and  the  reduced  normal 
momentum  equation  to  their  bleed  boundary  condition,  with  the  latter  case  yielding  a  better 
match  to  the  test  data  as  expected.  Rizetta^*  mentions  that  the  inclusion  of  the  normal  pressure 
gradient  equation  resulted  in  "smoother  behavior  of  the  flow  variables"  than  the  zero  normal 
pressure  gradient  assumption.  Benhachmi  et  al.^  also  used  the  reduced  normal  momentum 
equation  for  their  bleed  boundary  condition  study. 

Large  bleed/bypass  slots  typically  found  at  the  throat  of  external  and  mixed  compression 
inlets  have  been  modeled  directly  in  two  or  three  dimensions,  with  the  slot  opening  and  plenum 
chamber  as  part  of  the  grid  geometry.  The  flowfield  near  the  slot  opening  involves  imbedded 
shocks,  subsonic  pockets,  and  rapid  turning  into  the  plenum.  The  slot  flow  and  the  inlet  subsonic 
duct  flow  are  highly  coupled,  which  is  critical  for  inlet  stability.  It  is  not  likely  that  a  wall  grid 
with  a  mass  efflux  boundary  condition  spanning  the  slot  opening  will  accurately  replicate  this 
complex  physical  situation.  As  an  example,  Fujimoto,  Niwa,  and  Sawada**  have  done  a  realistic 
mixed  compression  inlet  calculation  including  ramp  slot  bleeds  and  a  throat  bypass  slot  with  the 
plenums  and  exhaust  ports  included.  Some  of  their  results  were  shown  in  Figure  2-11.  No  test 
data  is  available  but  this  is  an  excellent  illustration  of  the  complexity  involved  in  a  good  design 
application. 

4.4  Supersonic  Outflow  Boundaries 

Supersonic  outflow  boundaries  for  inlet  applications  are  those  which  involve  supersonic 
flow  leaving  the  computational  domain.  The  desired  condition  is  that  any  waves  (shocks,  Mach 
waves,  etc.)  pass  through  the  boundaries  without  reflection  in  accordance  with  the  direction  of 
the  flow  characteristics  at  that  boundary.  This  is  very  beneficial  because  these  boundaries  can 
be  placed  close  to  the  region  of  interest  without  influencing  the  solution  and  therefore  reduce  the 
domain  size.  Also,  mixed  compression  inlets  which  maintain  supersonic  flow  throughout  the 
internal  duct  (such  as  the  P8)  are  simulated  with  a  supersonic  outflow  boundary  at  the  duct  exit. 

The  nonreflective  property  is  accomplished  by  extrapolating  the  flow  variables  at  the 
boundary.  The  implementation  of  extrapolation  boundary  conditions  is  straightforward  by  two 
methods.  The  first  method  is  to  simply  specify  zero  gradients  in  the  conservative  variables  in 
a  direction  normal  to  the  boundary: 
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This  is  numerically  applied  with  a  two  or  three  point  backward  difference.  The  application 
references  use  different  terminology  to  describe  this  condition  without  giving  the  details, 
including  "nonreflective  boundary",  "extrapolation  boundary",  and  "zero  flux  gradient." 

The  second  method  is  a  characteristics  based  boundary  condition.  The  Riemann  invariants 
for  unsteady  wave  motion, 

AP+(u±a)Ap=0  4.6 


are  used  to  extrapolate  the  boundary  points  along  the  characteristic  directions.  This  condition  is 
based  on  the  physics  of  the  flow  but  more  complex  to  implement  numerically.  The  direction  of 
the  characteristics  and  their  intercept  points  on  the  boundary  must  be  determined  and  may  require 
an  iterative  process.  Terminology  for  this  condition  is  "Mach  wave  extrapolation"  or 
"characteristic  boundary." 

The  zero  flux  gradient  condition  has  been  the  most  widely  used  according  to  the 
application  survey,  probably  due  to  its  simplicity  and  satisfactory  results.  Both  methods  work 
well  as  long  as  the  boundary  is  a  reasonable  distance  away  from  high  flow  gradients  in  the  region 
of  interest.  "Reasonable"  means  the  region  of  interest  is  not  influenced  by  the  boundary,  which 
may  require  some  numerical  experimentation. 

4.5  Subsonic  Outflow  Boundaries 

A  boundary  where  subsonic  internal  duct  flow  exits  an  inlet  configuration  needs  special 
consideration.  The  exit  pressure  and  massflow  rate  are  critical  parameters  for  determining  the 
inlet  operating  point  and  performance.  Pressure  and  massflow  depend  upon  the  downstream 
device  which  is  throttling  the  inlet  (an  engine  for  air  vehicles,  a  massflow  plug  for  test  hardware, 
etc.).  Actually,  the  boundary  itself  is  a  computational  artifice  because  the  throttling  device  is  too 
complicated  to  simulate  in  a  CFD  solution  or  the  physical  situation  may  not  be  well  defined. 

High  pressure  subsonic  exit  conditions  impose  a  strong  upstream  influence  on  the  inlet 
shock  structure  which  affects  boundary  layer  interactions,  massflow  control,  and  inlet  stability. 
To  simulate  these  effects,  subsonic  exit  conditions  must  propagate  information  from  the  boundary 
to  the  interior  grid  points.  From  1-D  unsteady  wave  theory,  subsonic  flow  has  two  positive  slope 
characteristic  propagating  from  the  interior  towards  the  boundary  and  a  negative  slope 
characteristic  propagating  from  the  boundary  towards  the  interior.^^^®’^’  Thus  one  flow  parameter 
may  be  specified  on  the  boundary  and  the  others  are  extrapolated  from  the  interior. 

The  most  common  approach  is  to  specify  a  uniform  static  pressure  along  this 
boundary.^^’^®’^'*’'*^’'**  The  pressure  level  determines  the  inlet  operating  point  and  massflow  rate  in 
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the  steady-state  solution.  Alternatively,  a  certain  massflow  rate  (corrected  or  actual)  may  be 
desired  at  the  exit.  This  is  accomplished  by  integrating  the  boundary  after  each  iteration  to 
determine  the  actual  massflow  rate.  Then  the  pressure  level  is  adjusted  as  the  solution  develops 
to  rectify  any  difference  from  the  desired  massflow  rate,  but  the  effect  is  still  a  uniform  exit 
pressure. 

One  drawback  to  the  uniform  static  pressure  boundary  is  that  the  flow  entering  an  engine 
or  massflow  control  device  is  not  constrained  to  have  uniform  static  pressure.  This  can  be  a 
problem  if  the  diffuser  flow  path  has  relatively  high  curvature  near  the  exit.  Two  methods  have 
been  used  to  avoid  this  problem.  The  first  is  to  add  a  constant  area  straight  section  of  grid  to 
the  diffuser  exit  and  specify  the  uniform  pressure  level  at  the  exit  of  the  straight  section,  as 
shown  in  Figure  4-3.  The  length  of  the  section  may  be  on  the  order  of  a  few  duct  diameters  or 
duct  heights,  which  helps  remove  the  uniform  exit  pressure  constraint  from  the  point  of  interest. 
The  disadvantage  is  that  additional  grid  points  are  needed  for  the  solution.  The  other  approach 
is  to  add  a  "computational"  nozzle  to  the  diffuser  exit.  This  is  a  convergent-divergent  section 
added  to  the  diffuser  exit  which  controls  massflow  rate  by  the  size  of  the  throat.  The  nozzle  exit 
boundary  is  an  extrapolation  condition  as  discussed  in  Section  4.4  since  the  flow  is  reaccelerated 
to  a  supersonic  Mach  number.  This  is  a  realistic  simulation  because  it  throttles  the  inlet  in  a 
similar  way  as  an  engine  exhaust  nozzle  or  massflow  plug.  However,  it  can  also  be 
computationally  expensive.  Before  a  solution  is  obtained,  the  nozzle  throat  area  must  be  sized 
from  estimates  of  total  pressure  loss  through  the  inlet  and  an  estimated  nozzle  discharge 
coefficient.  Any  fine  tuning  of  the  throat  area  to  get  a  desired  massflow  rate  requires  both  a  new 
nozzle  grid  and  a  restart  of  the  solution.  Shigematsu  et  al.'*^  obtained  2-D  solutions  of  a  mixed 
compression  inlet  with  various  nozzle  exit-to-throat  area  ratios  as  shown  in  Figure  4-4.  In  this 
case  the  "computational"  nozzle  is  handled  by  their  flow  solver  such  that  only  the  throat  area 
needs  to  be  specified;  the  nozzle  grid  is  automatically  regenerated  internally  by  the  code.  The 
decreasing  nozzle  throat  area  results  in  enlarged  separated  regions  in  the  diffuser  as  the  increased 
back  pressure  feeds  forward.  Experience  suggests  that  similar  events  would  occur  in  a  wind 
tunnel  test. 
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Figure  4-3  Gridding  Techniques  for  Diffuser  Exit  Boundaries 


Figure  4-4 


Mixed  Compression  Inlet  with  Computational  Nozzle 
Exit  Boundary  (Ref.  42) 


5.0  GRID  REQUIREMENTS 


Grid  generation  appears  to  rely  heavily  upon  personal  experience  with  CFD  applications. 
There  are  few  guidelines  when  starting  a  new  application.  Often  some  adjustments  are  made 
after  an  initial  solution  is  obtained  to  improve  the  grid  quality  in  certain  regions.  Complexity  of 
the  grids  for  applications  reviewed  here  range  from  simple  single  block  rectangular  geometry  for 
shock  wave-boundary  layer  interactions  to  complex  multiblock  arrangements  for  complete  inlet 
hardware. 

The  importance  of  grid  orthogonality  at  the  boundaries  has  already  been  discussed  in 
Section  4.2.  Also,  smooth  variation  of  the  grid  spacing  is  desirable  to  integrate  efficiently  to  a 
converged  steady-state  solution.^^  Abrupt  changes  in  grid  spacing  leads  to  large  truncation  errors 
and  possible  departure.  Closely  spaced  grid  points  are  desirable  in  regions  with  high  gradients, 
e.g.  shocks  and  boundary  layers,  to  reduce  numerical  truncation  errors.  However,  computer 
memory  limitations  and  processor  requirements  can  restrict  the  total  number  of  grid  points.  Thus 
the  CFD  applier  must  trade  accuracy  against  computer  costs  when  generating  grids.  More  often 
than  not,  geometric  simplifications  are  made  to  reduce  grid  complexity  as  much  as  possible. 

A  common  example  of  geometric  simplification  is  the  assumption  of  perfectly  sharp  ramp 
and  cowl  leading  edges.  The  assumption  of  a  sharp  leading  edge  greatly  simplifies  the  gridding 
task  and  eliminates  an  otherwise  densely  packed  grid  around  a  blunt  leading  edge.  The  P8  inlet 
had  a  cowl  lip  radius  of  0.045  inch,  which  definitely  displaced  the  cowl  shock  forward  of  the  lip 
highlight  station.  Knight^’  simulated  this  displacement  with  a  sharp  cowl  lip  1.1  inches  forward 
of  the  actual  lip  highlight  such  that  the  cowl  shock  would  strike  the  ramp  shoulder  at  the  right 
location  (Figure  5-1).  However,  the  entropy  layer  created  by  the  blunt  lip  is  ignored  by  this 
technique. 

5.1  Shock  Wave  Resolution 

State-of-the-art  integration  schemes  are  currently  capable  of  capturing  shocks  with  three 
to  four  grid  points.  The  most  recent  scheme  development  has  been  focused  on  capturing  shocks 
without  producing  dispersive  errors  (post  shock  oscillations).  However,  many  "production"  codes 
being  applied  to  inlets  still  can  produce  significant  post  shock  oscillations  because  they  are  based 
on  less  advanced  integration  schemes.  If  the  grid  spacing  happens  to  be  large  where  a  shock  is 
located,  then  the  shock  will  be  smeared  over  a  non-physically  large  distance,  which  might  affect 
critical  shock  boundary  layer  interactions.  Successful  CFD  appliers  make  an  estimate  of  shock 
locations  and  construct  the  grid  accordingly  to  yield  crisp  shock  capturing  with  minimal 
overshoot.  Also,  adaptive  grid  techniques  have  been  used  to  pack  the  grid  near  shocks  as  the 
solution  develops  to  a  final  steady-state. 

5.2  Boundary  Layer  Resolution 

Turbulent  boundary  layer  resolution  also  demands  high  grid  density  near  solid  walls.  The 
parameters  used  to  describe  the  grid  spacing  at  the  wall  are  closely  related  to  turbulent  boundary 
layer  theory  and  the  dimensionless  inner  variables: 
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where  x^,  and  are  the  kinematic  viscosity,  shear  stress,  and  density  at  the  wall, 
respectively.  A  typical  incompressible  turbulent  boundary  layer  profile  in  these  coordinates  is 
shown  in  Figure  5-2  with  the  different  regions  identified.  The  location  of  the  first  grid  point 
from  the  wall  is  important  to  get  accurate  wall  shear  stress  and  heat  transfer.  This  location  is 
measured  in  terms  of  the  wall  normal  coordinate  y"^.  Since  y'^  depends  upon  the  local  wall  shear 
stress,  the  applier  will  not  know  the  exact  value  for  the  first  grid  point  prior  to  obtaining  the 
solution.  Estimates  from  flat  plate  profiles  are  used  to  construct  the  initial  grid,  with  adjustments 
made  as  the  solution  progresses.  The  final  quality  of  the  grid  is  often  reported  in  terms  of  all 
the  grid  points  next  to  solid  walls  being  less  than  a  certain  y*  threshold.  A  satisfactory  value  for 
y*  is  a  matter  of  debate,  and  depends  upon  the  assumptions  of  the  turbulence  model  being  used 
(Section  6.0).  Wall  shear  stress  and  heat  transfer  accuracy  improves  with  closer  spacing.  Also, 
it  is  generally  known  that  accurate  heat  transfer  prediction  demands  closer  spacing  than  accurate 
shear  stress  prediction.  But  extremely  small  spacing  can  severely  restrict  the  allowable 
integration  time  step  for  stability  which  results  in  excessive  iterations  to  convergence. 

From  the  reference  applications,  values  were  reported  from  y^=.07  to  as  much  as  y^=8, 
with  most  reported  values  in  the  range  y^=l-5.  A  general  observation  from  the  references 
reviewed  here  and  discussions  at  technical  conferences  and  symposia  is  that  y*=l  is  considered 
to  be  quite  good.  The  perception  is  that  the  applier  did  not  sacrifice  boundary  layer  resolution 
in  order  to  get  an  economic  solution.  Agarwal*’'*^  and  co-authors  suggest  that  y*  should  be  at 
most  one  for  accuracy  in  skin  friction  calculations.  This  conclusion  is  supported  by  numerical 
experiments  which  showed  that  skin  friction  becomes  independent  of  normal  grid  spacing  when 
y*  is  less  than  one. 

Some  turbulence  models  require  use  of  a  wall  function.^'*’^®’^”  The  wall  function  is  an 
analytical  expression  for  the  turbulent  boundary  layer  profile  between  the  wall  and  the  first  grid 
point  from  the  wall.  A  typical  expression  is  the  familiar  logarithmic  law-of-the-wall: 


u*  =  —  In  y*+B  ,  ic«.41  ,  B«5 
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5.2 


To  use  a  wall  function,  the  solution  is  integrated  on  the  interior  points  as  usual,  but  boundary 
conditions  are  applied  at  the  first  point  off  the  wall  rather  than  at  the  wall.  These  boundary 
conditions  are  related  to  the  wall  shear  stress  via  the  wall  function.  The  first  point  off  the  wall 
should  be  located  in  the  logarithmic  region,  such  that  much  larger  values  for  y^  can  be  used  with 
corresponding  increases  in  stable  time  integration  step  sizes.  For  the  applications  referenced 
(24,28,50),  wall  normal  grid  spacing  up  to  y'’=30  was  used.  Using  a  wall  function  assumes  that 
the  logarithmic  near  wall  profile  exists  throughout  the  flow.  This  is  not  likely  for  inlet  flows 
with  strong  shock  wave-boundary  layer  interactions,  and  therefore  the  wall  functions  have  not 
seen  much  application  in  this  area. 

The  number  of  grid  points  in  the  turbulent  boundary  layer  is  also  reported  in  application 
references  as  a  measure  of  refinement  at  solid  walls.  These  values  ranged  from  10  to  40  points 
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Figure  5-1 
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Figure 


Sharp  Lip  Extension  on  the  P8  for  Grid 
Simplification  (Knight,  Ref.  29) 


5-2  Turbulent  Boundary  Layer  Nomenclature  (Ref.  3) 
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in  the  boundary  layer,  tending  toward  the  higher  numbers.  Furthermore,  the  number  of  points 
in  the  viscous  sub-layer  was  also  reported  in  the  range  from  2  to  18.  Again,  the  applier  will  not 
know  beforehand  the  actual  boundary  layer  thickness,  so  estimates  must  be  used  initially. 

Some  appliers  carefully  investigate  the  effect  of  grid  point  density  on  the  simple  unit 
problems  before  attempting  a  more  complicated  configuration.  Such  problems  are  small  enough 
that  several  solutions  can  be  obtained  to  optimize  the  grid  spacing.  It  is  generally  true  from  the 
application  references  that  grid  refinement  in  the  boundary  layer  is  diminished  as  larger  and  more 
complicated  configurations  are  attempted.  This  is  evidence  that  computer  economics  is  a  strong 
driver  in  grid  generation. 
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6.0  TURBULENCE  MODELS  AND  TRANSITIONS  SIMULATIONS 


High  speed  inlet  flows  of  interest  are  mostly  turbulent.  The  Boussinesq  assumption 
discussed  in  Section  2.2  introduces  the  proportionality  coefficient,  /x,,  but  does  not  offer  any 
information  abouts  its  numerical  value.  The  eddy  viscosity  coefficient  is  not  a  property  of  the 
fluid  like  molecular  viscosity.  Rather,  it  is  a  property  of  the  flowfield  and  depends  upon  the 
flow  conditions  and  geometry.^  Turbulence  models  attempt  to  simulate  this  dependency  through 
analytical,  empirical,  and  heuristic  reasoning.  The  most  common  models  encountered  in  this 
review,  their  modifications,  and  successfulness  of  the  applications  are  discussed  below. 

6.1  Zero-Equation  Models 

The  zero-equation  models  are  referred  to  as  algebraic  models.  None  of  the  parameters 
in  the  model  are  additional  dependent  variables  to  be  solved  with  Equations  2. 1-2.4.  They  are 
merely  evaluated  with  algebraic  closure  expressions  that  depend  on  the  local  flow  conditions. 
They  are  relatively  fast  and  require  a  negligible  amount  of  additional  computer  operations. 
However,  they  are  also  strongly  empirical  and  therefore  not  universally  applicable  without 
modifications  for  specific  types  of  flows.  It  will  be  shown  that  the  existing  algebraic  models 
have  serious  deficiencies  for  high  speed  inlet  applications.  However,  they  will  continue  to  be 
employed  because  of  their  relative  simplicity  and  widespread  incorporation  into  existing  Navier- 
Stokes  solvers. 

The  Cebeci-Smith  and  the  Baldwin-Lomax  models  have  been  the  primary  algebraic 
models  used  in  high  speed  inlet  applications.  They  are  both  two-layer  models  with  different 
formulations  for  the  inner  and  outer  regions  of  the  boundary  layer: 

_  r  ^t.iaaeTf  Y^Yc  6.1 

^  t  l^t,  outer'  Y^Yc 

where  is  the  "crossover"  point  separating  the  implementation  of  the  two  formulation-  The 
outer  formulation  is  used  down  to  the  "crossover"  point  y^,  where  it  first  predicts  a  lowe.  eddy 
viscosity  than  the  inner  formulation.  Then  the  inner  formulation  is  invoked. 

For  both  models,  the  inner  formulation  is  based  on  Prandtl's  mixing  length  theory^: 

6.2 
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where  I  is  the  mixing  length,  represented  empirically  by 

1  =  xyD  (k  =  0.41)  6.3 


and  D  is  the  Van  Driest  damping  factor: 
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The  damping  is  needed  to  reduce  the  eddy  viscosity  in  the  viscous  sublayer  where  molecular 
viscosity  dominates.  The  derivative  term  in  Equation  6.2  is  the  magnitude  of  vorticity  in  a  2-D 
boundary  layer,  which  is  often  replaced  by  the  complete  vorticity  expression  in  3-D  problems: 

“  “  2 -D  boundary  layer 

6.5 

rx  -  (  du  dwy,  [  dv  9vv  ■  ,  n 


Up  to  this  point,  the  Cebeci-Smith  model  and  the  Baldwin-Lomax  model  are  identical. 
They  differ  in  the  formulation  for  the  outer  region.  For  the  Cebeci-Smith  model,  Clauser's 
formula  for  jets  and  wakes  is  used^^’^: 

1*  t,  outer  =  0-016  8p  1^5  *^icleb 
where  is  the  boundary  layer  edge  velocity. 


is  the  kinematic  displacement  thickness,  and 


1+5.5 

6* 


is  the  Klebanoff  correction  for  turbulent  boundary  layer  edge  intermittency.  The  drawback  with 
this  formulation  is  that  the  outer  length  scale  depends  upon  the  boundary  layer  thickness  and  edge 
velocity  which  are  not  known  a  priori  and  might  be  difficult  to  determine  for  shock  wave¬ 
boundary  layer  interactions.^  To  eliminate  this  problem,  the  Baldwin-Lx)max  model^’^  replaces 
Clauser's  wake  model  with 


t,  outer  0  •  0 16  8  p  CicjpJ^max^max'^JCleab 


Hit 

where  =  max[y|-^|D]  and  is  the  value  of  y  at  The  Van  Driest  damping 

factor  is  from  Equation  6.4  and  the  Klebanoff  intermittency  correction  is  now 
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The  Cebeci-Smith  model  and  the  Baldwin-Lomax  model  will  produce  similar  results  for  flat  plate 

incompressible  boundary  layers  when  Cg,=1.6  and  C^fci,=0-3. 

Note  that  the  above  developments  and  empirical  constants  are  based  on  incompressible 
turbulent  boundary  layers  in  the  absence  of  strong  pressure  gradients  or  wall  transpiration. 
Algebraic  models  have  difficulty  predicting  separation  and  reattachment  induced  by  adverse 
pressure  gradients  because  they  are  equilibrium  models.  This  assumes  that  the  rates  of  turbjlent 
energy  production  and  dissipation  are  equal  throughout  the  flowfield.  In  their  unmodified  term, 
they  cannot  simulate  lag  or  "history"  effects  through  shock  wave-boundary  layer  interactions. 
For  compressible  turbulent  interactions  with  bleed,  blowing,  etc.,  modifications  and  "tuning"  are 
necessary  to  get  reasonable  results.  These  modifications  are  discussed  below. 

A  number  of  modifications  for  compressibility  and  wall  suction  were  encountered  which 
involved  adjustments  to  the  von  Karman  constant  k  and  the  van  Driest  dainping  thickness  A  . 
Rizetta^^  simply  changed  these  constants  to  new  values  based  on  2-D  numerical  experiments  of 
his  Mach  5  inlet  application.  More  complicated  expressions  for  A*  were  used  by  Horstman  et 
aP  and  Abrahamson  and  Brower^  to  allow  variation  with  streamwise  pressure  gradient  and  wall 
suction  velocity.  Both  applications  were  of  shock  wave-boundary  layer  interactions.  Both  papers 
concluded  that  A*  had  no  impact  on  the  performance  of  the  Baldwin-Lomax  for  predicting 
reattachment  location.  This  was  attributed  to  the  limited  influence  of  the  van  Driest  damping 
factor  in  the  outer  layer,  where  is  orders  of  magnitude  greater  than  A*.  Visbal  and  Knight 
studied  the  Baldwin-Lomax  model  for  a  Mach  2.9  compression  ramp  shock-boundary  layer 
interaction.  They  indicated  that  and  in  the  outer  layer  formulation  were  functions  of  the 
freestream  Mach  number  (C,^=2.08  and  0.65  were  required  to  match  compressible  Mach 

3  flat  plate  turbulent  profiles).  A  number  of  researchers  have  adopted  their  suggestions  or 
investigated  new  expressions  to  allow  to  vary  with  edge  Mach  number  in  order  to  improve 
wall  pressure  and  skin  friction  predictions  for  strong  interactions. 

Visbal  and  BCnight  also  identified  a  fundamental  flaw  in  the  Baldwin-Lomax  model  when 
applied  to  separated  interactions.  In  a  non-separated  equilibrium  turbulent  boundary  layer,  the 
F  function  exhibits  a  single  distinct  maximum.  For  separated  regions,  however,  two  peaks 
appear.  As  the  flow  reattaches,  the  inner  peak  becomes  dominate  (Figure  6-1).  This  is  indicative 
of  the  incoming  turbulence  equilibrium  level  being  disrupted  by  the  sudden  pressure  gradient  and 
then  relaxing  to  a  new  equilibrium  state.  The  multiple  peak  profile  in  the  separated  region  makes 
ambiguous.  can  differ  by  an  order  of  magnitude  between  the  two  peaks;  selecting  the 
peak  closest  to  the’wall  yields  a  non-physical  drop  in  outer  eddy  viscosity.  Visbal  and  Knight 
corrected  this  situation  by  ensuring  that  the  outermost  peak  is  always  selected  for  F„^.  However, 
as  the  innermost  peak  in  the  F  function  becomes  dominant  during  reattachment,  it  remains  close 
to  the  wall  such  that  the  eddy  viscosity  becomes  non-physically  low  again.  Sakowski  et  al 
addressed  this  problem  by  fixing  upstream  of  the  interaction  and  determined  F„„  at  this 
constant  location  throughout  the  interaction.  Another  problem  was  identified  by  Visbal  and 
Knight  in  the  evaluation  of  the  van  Driest  damping  factor.  Since  D  depends  upon  wall  shear 
stress  through  y*,  the  damping  factor  will  become  very  small  when  wall  shear  goes  to  zero  at  the 
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points  of  separation  and  reattachment.  This  causes  another  non-physical  drop  in  outer  eddy 
viscosity.  To  alleviate  this  problem,  the  wall  shear  stress  is  replaced  by  the  local  shear  stress  at 
points  off  the  wall  in  evaluating  D.  This  correction  has  also  been  adopted  by  various 
researchers.^^^^’^^ 

An  equilibrium  relaxation  correction  to  algebraic  models  was  investigated  by  Shang  et 
al.^^”  and  Horstman  et  aP*  for  shock  wave-boundary  layer  interactions.  This  correction  attempts 
to  add  a  history  to  the  turbulence  level.  Experimental  evidence  suggests  that  Reynolds  stress  is 
convected  along  streamlines  at  constant  levels  through  a  sudden  adverse  pressure  gradient,  and 
then  exponentially  approaches  a  new  equilibrium  state.^  This  was  modeled  by  the  following 
expression: 

_  Ax  6  11 
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where  is  the  value  of  eddy  viscosity  upstream  of  the  interaction,  is  the  value  of  eddy 
viscosity  downstream  of  the  interaction,  \  is  the  relaxation  length  scale,  and  Ax  is  the  distance 
from  the  upstream  location.  This  correction  provides  remarkable  improvement  in  the  prediction 
of  the  extent  of  upstream  influence  and  the  overall  interaction  length  for  shock  wave-boundary 
layer  interactions.  The  value  of  \  is  probably  a  function  of  upstream  Mach  number  and  shock 
strength,  ranging  from  one  to  twenty  boundary  layer  thicknesses  in  the  references  discussed  here. 
A  major  drawback  of  this  correction  is  that  the  position  of  the  interaction  must  be  known 
beforehand,  which  limits  its  usefulness  for  inlet  design  applications. 

One  area  where  the  algebraic  models  consistently  fail  regardless  of  corrections  is  the 
reattachment  and  recovery  of  separated  regions.  The  skin  friction  is  consistently  predicted 
incorrectly.  Ong  and  KnighP  and  Visbal  and  KnighP  attribute  these  failures  to  a)  inability  of 
current  CFD  techniques  to  simulate  post-shock  turbulent  fluctuation  amplification  and  natural 
shock  unsteadiness,  and  b)  inability  of  the  algebraic  models  to  simulate  dual  length  scales  in  a 
reattaching  flow  (one  for  the  outer  wake,  another  for  the  newly  developing  boundary  layer). 
Based  on  these  conclusions,  one  cannot  expect  the  algebraic  models  discussed  here  to  provide 
accurate  flowfield  results  for  a  complex  high  speed  inlet.  The  algebraic  models  have  difficulty 
with  the  reattachment  region,  and  shock  wave-boundary  layer  interaction  characteristics  depend 
upon  the  condition  of  the  approaching  boundary  layer.  Thus,  mixed  compression  inlets  with 
multiple  shock  interactions  become  progressively  corrupted  as  the  solution  moves  downstream. 
This  is  evident  in  the  results  of  P8  solutions  and  Rizzeta's^^  Mach  5  inlet  solution.  Still,  inlet 
designers  continue  to  accept  the  limitations  of  the  algebraic  models  because  of  the  economical 
routine  application  they  afford. 

There  is  an  interesting  aspect  to  the  application  of  the  algebraic  models  for  inlets.  The 
eddy  viscosity  depends  upon  the  distance  from  the  wall.  But  inlets  may  have  more  than  one  wall 
in  close  proximity  to  the  point  of  interest.  Internal  passages  have  corners  and  opposing  walls 
with  interacting  turbulent  wake  regions.  It  is  possible  that  the  dominating  wake  region  for  a 
given  point  in  a  field  may  not  be  due  to  the  closest  wall.  A  number  of  different  approaches  have 
been  taken  to  account  for  the  influences  of  multiple  walls  on  the  eddy  viscosity  at  a  point  in  the 
field.  The  simplest  approach  was  to  simply  let  "y"  be  the  distance  to  the  closest  wall  and  ignore 
the  influence  of  other  walls. 
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Figure  6-1  Baldwin-Lomax  Model  F  Function  Peaks  in  a 

Separated  Region  (Ref.  23) 


Plane  of 
Symmetry 


Figure  6-2  Wall  Normal  Distances  for  Turbulence  Model 

Applications  in  Multi-wall  Configurations 
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35 

Other  approaches  involve  weighting  functions  based  on  distance  to  all  walls.  Rizetta 
used  the  following  formula  for  eddy  viscosity  at  a  point  in  the  symmetric  rectangular  duct  of  the 
Mach  5  inlet: 

_  6.12 

1/t1i+1/T]2+1/t13 

The  enumerated  eddy  viscosity  coefficients  are  obtained  independently  for  each  wall,  and  t]^,  tIj, 
and  rij  are  the  normal  distances  to  each  wall,  Figure  6-2.  Rizetta  admits  that  there  is  no  physical 
basis  for  this  formula,  but  it  gives  the  correct  form  away  from  the  comers  and  is  continuous.  A 
variation  on  this  formula  is  given  by  Ramakrishnan  and  Goldberg^’: 


= 
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6.13 


This  formula  places  more  weight  on  the  closest  wall,  but  again  there  was  no  physical  basis  for 
this  discussed  in  the  reference.  Knight^®’'‘°'^^  and  Shang  et  al'^^  have  used  a  more  general  integral 
form  due  to  Buleev: 


where  s  is  the  distance  from  the  point  of  interest  to  a  point  on  the  wall  boundary  defined  by  the 
angle  (j)  (Figure  6-3).  This  integral  results  in  algebraic  expressions  for  orthogonal  corners  and 
rectangular  ducts.  The  length  I  is  then  used  in  the  algebraic  eddy  viscosity  model  in  place  of  "y". 
The  Buleev  length  scale  is  apparently  based  upon  theory  and  experimental  observation  of 
turbulent  corner  flows. 

Still  another  method  was  used  by  Ng  et  al.^^  for  a  2-D  analysis  of  the  P8  inlet.  In  this 
situation  the  lower  wall  will  have  a  much  thicker  wake  region  than  the  upper  wall,  and  thus  the 
lower  wall  could  have  more  influence  on  the  eddy  viscosity  at  points  which  are  physically  closer 
to  the  upper  wall.  Their  method  is  a  two  step  process  which  first  determines  two  independent 
values  of  eddy  viscosity  at  a  point  for  both  the  upper  and  lower  walls.  Then  these  values  are 
compared  and  the  greater  value  is  selected  to  advance  the  solution. 
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6.2  Two-Equation  Models 


The  two-equation  models  are  more  complex  than  the  algebraic  models,  but  they  still 
invoke  the  Boussinesq  assumption  and  relate  turbulent  Reynolds  stresses  to  the  flow  gradients 
via  the  eddy  viscosity  coefficient.  The  complexity  arises  from  a  functional  relationship  between 
the  eddy  viscosity  and  two  additional  dependent  conservation  variables  that  describe  the  velocity 
and  length  scales  of  the  turbulence.  The  two  additional  dependent  variables  have  their  own 
transport  differential  equations  with  terms  for  convection,  diffusion,  dissipation,  and  production. 
In  this  way  the  eddy  viscosity  is  related  to  the  flowfield  gradients  with  minimal  empirical 
presumptions  about  the  shape  of  the  boundary  layer  profiles.^  The  major  advantage  advertised 
for  the  two-equation  models  is  their  ability  to  include  upstream  influence,  or  the  so-called 
"history"  of  turbulent  stresses.^’ Recall  that  the  failure  of  the  algebraic  models  for  shock 
wave-boundary  layer  interactions  was  attributed  in  part  to  the  assumption  that  the  turbulence  was 
in  equilibrium.  The  two-equation  models  do  not  make  this  assumption:  production  and 
dissipation  rates  of  turbulent  energy  may  be  unequal.  Therefore  one  might  expect  the  two- 
equation  models  to  perform  better  for  such  interactions. 

Most  of  the  relevant  applications  of  two-equation  models  to  high  speed  inlet  flows  were 
shock  wave-boundary  layer  interaction  unit  problems.  The  only  application  found  for  an  actual 
inlet  with  test  data  was  the  P8  inlet.  This  attests  to  the  difficulties  of  applying  the  two-equation 
models  in  a  full  Navier-Stokes  analysis  of  a  complete  inlet  geometry.  The  two  additional 
differential  equations  add  considerable  computer  memory  and  processing  requirements  to  state-of- 
the-art  implicit  flow  solvers.  With  the  need  for  very  fine  grid  resolution  near  solid  walls,  most 
complete  inlet  analysis  problems  quickly  exceed  available  computer  resources.  Also,  the 
turbulence  transport  equations  are  usually  very  "stiff"  (i.e.,  their  solution  relaxes  slowly  from 
initial  conditions  to  a  time  asymptotic  solution),  which  requires  a  large  number  of  iterations  to 
arrive  at  a  converged  steady-state  solution. 

One  method  of  reducing  the  required  computer  resources  (memory  and  operations)  is  to 
decouple  the  two  turbulence  transport  equations  from  the  Navier-Stokes  equations  by  time 
lagging.  For  each  global  iteration  in  time.  Equations  2.1-2.4  are  first  solved  with  values  for  the 
two  turbulence  conservation  variables  (and  thus;^,)  from  the  previous  iteration.  Then,  the  two 
turbulence  transport  equations  are  solved  by  themselves  using  the  results  of  the  first  step  to 
update  the  turbulence  conservation  variables  for  the  next  global  iteration. 

The  two-equation  models  encountered  in  the  literature  search  involve  the  turbulent 
conservation  variables  k,  £,  and  co.  These  variables  generally  refer  to  turbulent  kinetic  energy, 
turbulent  dissipation  rate,  and  specific  turbulent  dissipation  rate,  respectively.  Relationships 
between  them  and  the  turbulent  length  scale  /  are  as  follows:^’^'*’^^ 
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The  eddy  viscosity  is  found  from  the  turbulent  conservation  variables  via  the  following  relations: 


V-t 
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General  forms  for  the  turbulent  transport  equations  are: 


=  convec ti on  +  di ffusi on  +  produc ti on  +  di ssipa tion  6.17 


The  terms  for  diffusion,  production,  and  dissipation  involve  various  combinations  of  flow 
gradients  and  algebraic  constants  which  are  specific  to  each  model.  The  combination  of  the  k 
transport  equation  with  either  of  the  others  yields  the  full  two-equation  model  description. 

A  particular  model  is  usually  associated  with  the  name(s)  of  its  originator.  The  following 
table  lists  the  models  which  were  encountered  in  the  literature  search: 


Two-Equation  Models  Encountered 
in  the  Literature  Search 


Name 

Type 

Notes 

Chien 

k-E 

low  Reynolds  number 

Jones-Launder 

k-E 

low  Reynolds  number 

Speziale 

k-E 

low  Reynolds  number 

Launder-Spaulding 

k-E 

high  Reynolds  number 

Wilcox-Rubesin 

k-o) 

In  the  low  Reynolds  number  form,  the  two  transport  equations  can  be  integrated  right  down  to 
the  wall,  with  for  the  first  grid  point  adjacent  to  the  wall  on  the  order  of  High 

Reynold's  number  formulations  do  not  model  effects  near  the  wall  adequately  and  require  the  use 
of  a  wall  function.®’^"*’^®  Typical  y*  values  for  the  first  grid  point  are  on  the  order  of  thirty.  Each 
model  has  different  assumptions  and  parameters  to  derive  the  terms  in  the  transport  equations. 
Little  detail  is  provided  in  the  application  references  about  these  terms,  and  the  authors  refer  the 
reader  to  other  references  for  more  information.  Other  modifications  have  been  included  for 
compressibility,^®’^®  thin  layer  approximation,®  and  damping  of  the  resultant  eddy  viscosity  in  the 
viscous  sublayer  where  molecular  viscosity  dominates.^''  Some  of  the  references  mention  that  the 
two-equation  models  are  prone  to  numerical  instabilities,  and  parameters  were  added  to  control 
them.2^’2®’^® 
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Additional  boundary  conditions  are  needed  since  two  new  equations  have  been  added  to 
the  system.  At  the  wall,  the  turbulent  kinetic  energy  is  assumed  to  be  zero,  and  a  zero  normal 
gradient  is  assumed  for  the  dissipation  variable.  This  is  consistent  with  the  laminar  sublayer  next 
to  the  wall  in  which  the  turbulence  is  damped  out.  At  an  exit  boundary,  both  variables  are 
extrapolated  from  the  interior.  At  an  entrance  boundary,  the  variables  are  held  fixed.  Agarwal, 
Deese,  and  Peters®  used  empirical  prescriptions  for  a  shock-boundary  layer  interaction  unit 
problem.  The  velocity  scale  is  related  to  freestream  Mach  number,  and  the  length  scale  is  related 
to  entrance  boundary  layer  thickness.  Lowrie^®  used  similar  empirical  prescriptions  for  a 
compression  ramp  problem,  but  cautioned  that  the  solution  results  were  sensitive  to  freestream 
length  scale.  Lowrie  further  suggests  that  it  is  appropriate  to  study  the  affect  of  freestream  length 
scale  for  new  applications. 

Many  of  the  applications  of  two-equation  turbulence  models  were  directed  toward  study 
of  model  performance  and  comparison  to  other  models.  The  collective  results  of  these  studies 
are  not  conclusive  about  the  superiority  of  one  model  over  another  or  two-equation  models  over 
algebraic  models.  Some  examples  are  cited  below. 

Agarwal,  Deese,  and  Peters®  compared  several  turbulence  models  for  various  shock¬ 
boundary  layer  interactions,  including  separated  and  non-separated  cases.  The  k-£  model  was 
compared  to  the  Baldwin-Lomax  model  for  a  separated  normal  shock  interaction  at  Mach  1.4. 
The  wall  pressure  distributions  predicted  by  these  two  models  were  identical.  Discrepancies  to 
the  test  data  show  that  the  upstream  influence  is  too  slight,  and  the  downstream  pressure  rise  is 
too  strong.  Both  models  predict  the  experimental  separation  and  reattachment  points  well  for  this 
case.  The  skin  friction  distribution  is  improved  by  the  k-s  model  downstream  of  reattachment, 
however. 

Carrol  and  Dutton^'^  studied  the  differences  between  the  Wilcox-Rubesin  model  with  a 
wall  function  and  the  Baldwin-Lomax  model  for  a  2-D  normal  shock  train  at  Mach  1.61.  The 
two-equation  model  was  superior  to  the  algebraic  model,  but  still  not  totally  in  agreement  with 
the  test  data.  The  authors  indicate  that  the  predicted  shock  train  for  both  models  had  a  tendency 
to  move  farther  downstream  in  the  channel  than  the  test  shows.  The  authors  arbitrarily  raised 
the  channel  exit  pressure  to  position  the  initial  shock  at  the  experimentally  observed  location. 
The  authors  suggest  that  boundary  layer  growth  on  the  end  walls  of  the  test  channel  (aspect  ratio 
of  2.4)  contributed  total  pressure  losses  which  resulted  in  a  lower  exit  pressure  in  the  test.  The 
authors  also  indicate  that  although  the  basic  flow  structure  was  not  affected  by  grid  density, 
shock  train  position  was  dependent  on  grid  refinement  in  the  transverse  direction. 

Kapoor,  Anderson,  and  Shaw^  applied  several  algebraic  models  and  all  the  k-E  models 
listed  in  Table  1  except  Jones-Launder,  to  the  P8  inlet.  The  calculations  were  fully  turbulent; 
boundary  layer  transition  was  not  simulated.  The  two  low  Reynolds  number  formulations  (Chien, 
Speziale)  matched  the  experimental  wall  pressures  through  the  cowl  shock  interaction  region  quite 
well.  One  of  the  algebraic  models  was  also  shown  to  provide  good  wall  pressure  distributions. 
The  high  Reynolds  number  formulation  (Launder-Spaulding)  over  predicted  the  peak  pressure  rise 
produced  by  the  interaction.  All  the  models  exhibited  too  much  upstream  influence  ahead  of  the 
cowl  shock  interaction.  The  authors  attribute  this  to  the  fact  that  the  calculations  were  fully 
turbulent.  This  yields  a  thicker  boundary  layer  than  the  experimental  conditions  with  transition, 
which  would  tend  to  move  the  Interaction  region  upstream. 
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With  regards  to  pitot  pressure  profiles  at  the  various  stations  in  the  P8,  the  k-s  models  are 
generally  closer  to  the  experimentally  observed  boundary  layer  profiles  than  the  algebraic  models. 
The  authors  concluded  that  it  is  acceptable  to  use  an  algebraic  model  to  predict  wall  pressures 
for  inlets  with  weak  interactions,  like  the  P8,  while  low  Reynolds  number  k-s  models  are  favored 
to  resolve  boundary  layer  details.  However,  the  review  of  the  Baldwin-Lomax  model  in  Section 
6.1  indicates  that  this  is  not  generally  true  for  strong  interactions. 

Lowrie  extended  the  low  Reynolds  number  Jones-Launder  model  to  a  multi-zone  grid 
approach  which  is  often  used  for  complex  configurations  to  simplify  grid  generation.  The 
turbulence  model  and  grid  zone  coupling  strategy  were  tested  against  a  24©  compression  ramp 
unit  problem  at  Mach  2.85,  which  produces  a  separated  interaction.  The  k-e  model  did  not 
accurately  predict  the  extent  of  the  separated  region,  regardless  of  whether  a  multi-zone  grid  was 
used  or  not.  The  author  attributes  this  failure  to  the  isotropic  assumption  of  the  eddy  viscosity 
concept  without  additional  analysis  or  discussion. 

6.3  Boundary  Layer  Transition  Models 

It  was  mentioned  earlier  that  inlet  flows  are  primarily  turbulent.  There  are  situations, 
however,  where  boundary  layer  transition  occurs  and  has  an  affect  on  inlet  performance  and  more 
importantly  inlet  surface  heat  transfer.  This  situation  is  common  in  wind  tunnel  tests  where  the 
Reynolds  number  is  low  compared  to  typical  flight  values.  This  was  the  case  for  the  P8  inlet, 
where  it  was  experimentally  observed  that  transition  occurred  at  about  40%  of  the  length  of  the 
ramp.  It  is  obviously  important  to  consider  transition  when  comparing,  validating,  verifying,  etc., 
CFD  methods  against  test  data. 

Methods  encountered  for  simulating  transition  in  the  inlet  applications  are  fairly  crude 
empirical  approximations  based  on  flat  plate  boundary  layer  experiments.  Most  of  the  methods 
require  a  specified  location  of  transition  onset  and/or  distance  of  the  transition  zone.  This  is  an 
extremely  limiting  requirement,  because  transition  details  are  seldom  known  or  cannot  be 
predicted  with  current  analysis  tools  for  complex  3-D  configurations. 

A  method  originally  developed  by  Dhawan  and  Narasimha  in  1958  was  outlined  by 
Shang,  Hankey,  and  Petty^^  and  implemented  on  a  supersonic  compressing  corner  at  Mach  3.  A 
transition  weighting  function  is  given  by 

r(^  =l-e3q3(-.4125^) 
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The  beginning  and  end  of  the  transition  region  are  specified  by  x^^^  and  x^^j,  respectively.  A  is 
a  measure  of  the  extent  of  the  region:  experimental  evidence  suggests  that  the  Reynolds  number 
at  the  end  of  transition  is  about  twice  that  at  the  start.  Further  details  of  the  implementation  were 
not  disclosed,  but  it  is  surmised  that  the  weighting  function  F  determines  the  fraction  of  eddy 
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viscosity  that  is  actually  added  into  the  governing  equations: 

total  “  l^Iaai  ^  l^inodel 


6.20 


In  this  implementation,  the  eddy  viscosity  from  the  turbulence  model,  must  be  determined 
throughout  the  flow  domain.  However,  only  a  fraction  is  added  to  the  laminar  molecular 
viscosity  per  the  exponentially  growing  function  F.  Note  that  this  method  does  not  actually 
predict  the  onset  of  transition;  it  only  describes  how  the  viscosity  varies  through  the  transition 
region.  Shang  et  al,  did  not  study  the  performance  of  this  model  for  their  compressing  corner 
problem. 

The  method  of  Dhawan  and  Narasimha  was  also  used  by  Knight  to  calculate  the  P8  inlet 
in  1977.^’  As  mentioned  above,  the  transition  location  of  this  inlet  was  experimentally  observed. 
The  detailed  affects  of  the  transition  model  were  not  specifically  addressed,  but  the  general 
results  in  terms  of  wall  and  pitot  pressure  comparisons  are  among  the  best  of  any  applications 
of  the  P8  uncovered  in  the  literature  search.  This  is  interesting  because  Knight  used  the  Cebeci- 
Smith  algebraic  turbulence  model.  Other  applications  of  the  P8  with  fully  turbulent  flow  using 
algebraic  models  produced  cowl  shock  induced  separation  bubbles  which  were  not  observed  in 
the  test.  This  suggests  that  the  turbulence  models  are  not  completely  at  fault,  and  that  the 
condition  of  the  boundary  layer  approaching  a  shock-boundary  layer  interaction  must  also  be 
simulated  correctly. 

Another  model  of  transition  was  used  on  the  P8  by  Yaghmaee  and  Roberts.^"*  In  this  case, 
the  model  is  capable  of  predicting  the  location  of  transition  onset  based  on  local  conditions: 

=  constant  6.20 

Me 


where  Reg  and  are  the  local  boundary  layer  momentum  thickness  Reynolds  number  and  edge 
Mach  number,  respectively.  The  constant  is  an  empirically  determined  threshold.  When  the  ratio 
on  the  left  exceeds  the  threshold,  the  turbulent  eddy  viscosity  is  ramped  up  linearly.  Another 
constant  is  specified  to  signal  the  end  of  transition,  at  which  point  the  full  value  of  eddy  vi^osity 
is  included.  This  method  is  particularly  useful  for  flows  experiencing  pressure  gradients  since 
it  is  based  on  local  conditions  and  has  been  recently  popular  for  engineering  studies  of  hypersonic 
vehicles.  Yaghmaee  and  Roberts  varied  the  onset  threshold  to  match  the  transition  point  for  the 
P8.  Again,  they  did  not  specifically  examine  the  affects  of  transition.  The  pitot  pressure  profile 
upstream  of  the  cowl  shock  interaction  compares  well  with  the  test  data,  indicating  that  the 
boundary  layer  approaching  the  interaction  was  simulated  correctly. 
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7.0  COMMENTS 


Based  on  the  discussion  of  Section  2.0,  it  can  be  concluded  that  the  state-of-the-art  in 
modeling  of  the  physics  of  high  speed  fluid  flow  is  sufficient  to  yield  an  accurate  simulation  of 
typical  inlet  configurations,  with  a  rather  restrictive  caveat  on  modeling  of  the  turbulent  boundary 
layer  interactions.  Commercial  CFD  codes  are  now  available  which  include  options  for  solving 
full  Navier-Stokes,  PNS,  and  Euler.  These  codes  include  perfect  gas,  equilibrium  high 
temperature  gas,  and  even  finite  rate  chemistry  models  with  broad  ranges  of  chemical  reactions 
from  which  to  select. 

The  advanced  numerical  integration  schemes  mentioned  in  Section  3.0  are  also  sufficient 
to  accurately  simulate  high  speed  inlet  flows.  Shock  capturing  capability  with  the  advanced 
schemes  is  not  a  concern,  although  some  experience  is  generally  necessary  in  order  to  execute 
these  advanced  methods  successfully.  Commercial  code  developers  often  claim  there  products 
are  "user  friendly,"  but  this  always  implies  that  new  problems  are  easy  to  set  up  and  initiate.  It 
does  not  mean  that  users  who  are  unfamiliar  with  aerodynamic  theory  and  numerical  methods  can 
apply  them  successfully. 

Grid  generation  tools  have  also  been  developed  by  government  and  industry  which  make 
this  tedious  task  more  automated  but  not  independent  of  the  skill  and  knowledge  of  the  user. 
Current  codes  and  grid  generation  tools  are  capable  of  multiblock  grid  solutions  which  are 
applicable  to  practically  any  supersonic  inlet  geometry.  The  complexity  of  an  inlet  application 
is  limited  by  the  size,  cost,  and  availability  of  computer  resources,  not  the  software.  New 
techniques  for  unstructured  grids  are  making  the  grid  generation  task  less  demanding.  These  are 
currently  being  applied  and  validated,  but  no  published  results  were  available  at  the  time  of  the 
literature  survey. 

Current  techniques  for  boundary  conditions  discussed  in  Section  4.0  cannot  be  considered 
mature  in  all  aspects.  The  typical  specifications  for  small  unit  problems  are  clearly  sufficient 
(freestream,  extrapolation,  solid  walls).  Issues  such  as  wall  pressure  normal  gradient  and 
temperature  conditions  have  either  minor  impact  or  their  influence  on  the  results  are  understood. 
However,  great  simplications  are  necessary  to  implement  numerical  boundaries  on  realistic  inlet 
configurations  in  order  to  limit  the  problem  domain  to  a  tractable  size.  Specifically,  the  bleed 
boundaries  and  diffuser  exit  boundaries  are  ncessarily  simplified,  with  questionable  validity  in 
some  cases.  New  techniques  for  treating  bleed  boundaries  appear  promising.  The  author  is 
aware  of  ongoing  research  in  this  area,  both  experimental  and  numerical.  Techniques  for 
improving  diffuser  exit  boundaries  are  also  under  development  which  simulate  the  flow  control 
characteristics  of  high  speed  engines. 

Turbulence  models  are  definitely  the  pacing  technology  issue  for  any  CFD  application 
which  involves  strong  shock  wave-boundary  layer  interactions.  This  issue  does  not  need  to  be 
belabored,  for  there  is  exhaustive  research  in  this  area.  In  the  mean  time,  continuous  application 
and  comparison  with  experiment  is  the  standard  operating  procedure. 
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